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INTRODUCTION 


The  flowing  medium  in  a  gun  tube  typically  is  a  mixture  of  a 
compressible  gas  with  burning  solid  propellant  grains.  Details  of  the  flow 
are  important  for  weapons  development,  but  only  bulk  properties  can  be 
routinely  measured,  such  as  the  trajectory  of  the  projectile,  the  pressure 
history  at  a  fixed  station,  the  heating  of  the  gun  tube,  etc.  Therefore,  a 
need  exists  for  a  detailed  mathematical  model  of  interior  ballistics  two- 
phase  flows. 

A  complete  mathematical  description  of  the  flow  could  provide  the 
motion  and  combustion  history  of  each  propellant  grain  and  of  the  gas  flow 
between  the  grains.  The  corresponding  local  governing  equations  are  easily 
established,  but  they  cannot  be  solved  numerically  because  of  the  great 
number  of  grid  points  needed  to  describe  a  flow  with  many  moving 
interfaces.  The  computational  work  can  be  reduced  only  be  sacrificing  the 
detailed  description  of  the  flow.  To  that  end  one  considers  mean  values  of 
the  two-phase  flow  that  are  derived  from  the  local  properties  of  the  gas  and 
grains.  The  governing  equations  for  these  average  properties  are 

established  by  averaging  the  local  governing  equations. 

This  report  presents  a  complete  and  consistent  mathematical  model  of 
three-dimensional,  transient  interior  ballistics  (gas-solid)  phenomena  in 
which  the  total  effects  of  the  gas  phase  viscosity,  turbulence,  and  heat 
conduction  on  the  average  variables  are  included.  In  contrast,  most 
existing  models  neglect  viscous  and  heat  conduction  effects,  and,  thus,  can 
characterize  only  the  wave  propagation  in  a  two-phase  flow.  The  theory  of 
the  model  is  complete  and  consistent  in  that  all  the  averaged  variables, 
equations,  initial  and  boundary  conditions,  regions  of  definition  of  the 
variables  and  correlations  are  precisely  defined  and  derived  using  the  same 
averaging.  The  need  for  such  an  approach  is  due  to  the  complexity  of  the 
multiphase,  multidimensional  viscous  flow  field  and  a  lack  of  detailed 
experimental  data.  Uhder  such  conditions,  models  formulated  on 

phenomenological  arguments  are  often  unreliable.  Also  a  phenomenological 
derivation  seldom  provides  precise  error  bounds.  A  theoretically  derived 
model  permits  one  to  investigate  with  more  confidence  ballistic  processes 
that  cannot  be  observed  in  detail  because  error  bounds  are  precisely 
formulated  and  can  be  tested.  Furthermore,  a  careful  mathematical 
derivation  of  the  model  can  reveal  restrictions  on  the  model  itself.  The 
presented  mathematical  model  possesses  the  following  features:  (1)  The 

averaging  process  insures  a  sufficient  differentiability  of  the  average 
variables  so  that  the  governing  partial  differential  equations  are 
defined.  (2)  Appropriate  averages  are  used  for  quantities  that  are  defined 
over  volume  and  for  those  that  are  defined  over  a  surface.  (3)  The  regions 
of  definition  of  the  average  variables  are  given.  (4)  The  necessary 
auxiliary  conditions  to  the  governing  equations,  e.g.  ,  initial  conditions 
and  boundary  conditions,  are  consistent  with  the  averaging  process  used  to 
derive  the  governing  equations.  (5)  Terms  that  are  modeled  by  correlations 
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possess  simple  physical  interpretations.  Estimates  are  given  for  the 
difference  between  the  theoretical  definitions  of  the  correlations  and  the 
expressions  actually  used.  (6)  Whenever  the  contribution  of  a  term  is 

neglected  in  an  equation,  a  corresponding  error  term  is  established.  (7) 
Because  the  equations  are  to  be  solved  numerically,  attention  is  given  to  an 

appropriate  form  for  numerical  solution.  (8)  The  model  represents  a  two- 

phase  (gas-solid)  flow  in  which  the  solid  ignites  and  burns,  and  it  also 
simulates  other  phenomena  which  occur  in  a  viscous,  heat  conducting  interior 
ballistic  flow. 

Previous  work  on  two-phase  equations  for  interior  ballistics  has  been 
done  by  Gough,  *  Kuo  et  al«^  Fisher  and  Trippe^  and  Krier  et  al.^  The 
primary  purpose  of  these  works  was  the  investigation  of  the  wave  propagation 
within  the  gun  tube  during  the  early  phases  of  the  interior  ballistic 

phenomenon.  Gough’s  equations  were  later  augmented  to  include  gas-phase 
viscosity  and  heat  conduction,  and  used  in  a  computer  program  developed  by 
Gibeling  et  al.  ^  Our  equations  are  different  because  we  have  used  a 
different  averaging  process,  chosen  a  different  set  of  dependent  variables, 
and  changed  some  correlation  models  that  provide  experimental  input  to  the 
theory.  Furthermore,  our  approach  differs  from  the  ones  mentioned  above 
because  it  is  based  solely  on  a  consistent  mathematical  theory. 

The  averages  in  this  report  are  computed  by  weighted  averaging  over  a 
finite  volume.  Gough*  used  instead  a  weighted  averaging  over  an  infinite 
space-time  domain  with  an  unspecified  weight  function.  The  rationale  of  our 
choice  is  based  on  the  observation  that  any  averaging  smooths  out  local 
details.  In  order  not  to  lose  too  many  details,  one  should,  therefore,  use 


P.5.  Gough,  " The  How  of  a  Compressible  Gas  Through  an  Aggregate  of  Mobile, 
Reacting  Particles ,ff  Ph.D .  Thesis,  Department  of  Mechanical  Engineering, 
McGill  University,  Montreal,  1974 . 

2 

K.K.  Kuo ,  J.H.  Koo ,  T.R.  Davis ,  and  G.R .  Coates ,  " Transient  Combustion  in 
Mobile,  Gas-Permeable  Propellants,"  Acta.  Astron. .  Vol.  3,  No.  7-8,  pp. 
5  74-591,  1976. 

3 

E.B.  Fisher  and  A.P.  Trippe,  "Mathematical  Model  of  Center  Core  Ignition  in 
the  175mm  Gun,"  Calspan  Report  VQ-516S-D-2,  1974. 

4N.  Krier,  W.F.  van  Hassell,  S.  Rafan,  and  J.  Vershau,  ’Model  of 
Flamespreading  and  Combustion  Through  Packed  Beds  of  Propellant  Grains, " 
University  of  Illinois  at  Urbana-Champaign  Report,  TR-AAE-74-1,  1974. 

^H.J.  Gibeling,  R.C.  Buggeln,  and  H.  McDonald,  "Development  of  a  Ttvo- 
Dimensional  Implicit  Interior  Ballistics  Code,  "  USA  ARDC  AMCCOM/Ballistic 
Research  Laboratory  Contractor  Report,  ARBRL-CR-00411,  APG,  MD,  January 
1980,  AD  No.  AD  387  458. 
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the  smallest  averaging  domain  that  is  compatible  with  the  requirements  of 
the  problem  at  hand.  One  requirement  of  the  averages  is  that  they  should  be 
differentiable  as  many  times  as  the  ensuing  governing  equations  indicate. 

It  has  been  shown  by  Delhaye  and  Achard6  that  line  or  surface  averages  of  a 
gas/particle  mixture  do  not  possess  the  required  differentiability 
properties.  Therefore,  the  smallest  domain  for  averaging  is  a  three- 
dimensional  volume.  Time  averaging  is  not  needed  to  insure 

differentiability,  if  the  weight  function  for  space  averaging  is  chosen 
properly  (see  Section  2.2).  If  one,  nevertheless,  chooses  to  time  average, 
then  the  time  average  interval  would  have  to  be  very  small  because  we  are 
interested  in  an  accurate  characterization  of  a  rapidly  changing  flow  field. 

The  size  of  the  averaging  volume  is  important.  The  use  of  an  infinite 
volume  for  averaging  is  not  appropriate  in  confined  flows  because  for  such  a 
volume  the  sura  of  the  volume  fractions  of  the  two  phases  is  not  equal  to 
one.  This  creates  problems  for  the  formulation  of  the  governing  equations 
and  the  boundary  conditions,  and  for  the  interpretation  of  the  results.  The 
problem  with  the  formulation  of  the  equations  is  eliminated  by  using  an 
appropriate  finite  volume  average,  while  the  others  become  more  easily 
tractable.  We  discuss  the  problems  in  Sections  4.4  and  4.6.  If  the  weight 
function  in  any  infinite  volume  average  is  zero  outside  some  finite  distance 
from  the  point  at  which  the  average  is  taken,  then  the  resulting  average  is 
obviously  equivalent  to  a  finite  volume  average.  If  the  value  of  the  weight 
function  is  zero  outside  some  distance  which  depends  on  the  location  of  the 
point  at  which  the  average  is  taken,  then  the  resulting  average  is 
equivalent  to  a  variable  finite  volume  average.  In  this  type  of  average 
additional  terms  in  the  partial  differential  equations  for  the  average 
quantities  appear  that  represent  the  effects  of  the  change  of  the  averaging 
volume  in  time  and  space.  This  complication  is  avoided  in  the  present 
report  by  restricting  the  attention  to  a  fixed  finite  volume  average  with  a 
fixed  weight  function. 

The  average  equations  which  are  derived  in  Section  3  include  the 
effects  of  gas  viscosity  and  of  turbulence.  Furthermore,  the  choice  of 
equations  for  averaging  and  the  choice  of  dependent  variables  has  a  bearing 
on  the  numerical  solution  of  the  equations.  We  have  chosen  a  set  of 
variables  that  eliminates  some  possible  numerical  singularities,  enhances 
the  accuracy  of  numerical  differentiation,  and  separates  important  physical 
processes  for  easier  modeling.  The  choice  of  variables  is  discussed  in 
Section  4.2.  We  also  have  chosen  the  internal  energy  equation  for  averaging 


J.M.  Delhaye  and  J.L.  Aehard,  "On  the  Use  of  Averaging  Operators  in  Two- 
Phase  Modeling f  in  Thermal  and  Hydraulic  Aspects  of  Nuclear.  Ueoatnn  Safety., 
Vol.  1:  Light  yqfrer  fleagtore ,  O.C.  Jones  and  S.G.  Bankoff ,  eds .,  pp.  289- 
332,  ASME,  New  York,  1977. 
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instead  of  the  commonly  used  total  energy  equation.  The  reasons  for  this 
choice  are  that  it  produces  a  clear  separation  of  physical  effects  and  a 

more  lucid  modeling  of  two- phase  phenomena.  They  are  discussed  in  Sections 
3.2.3  and  4.7.3,  respectively.  As  a  result  of  the  consideration  of  viscous 
effects  and  the  choice  of  equations  and  variables,  our  governing  equations 
differ  from  those  derived  by  Gough.  Each  set  of  equations  has  different 
appro  cimation  errors  and  some  of  the  required  models  of  experimental 
correlations  are  different. 

The  experimental  correlations  in  interior  ballistics  are  characterized 
by  a  scarcity  of  data.  This  is  one  reason  why  corresponding  mathematical 
models  have  not  been  firmly  established.  In  Section  4.7  we  list  a  set  of 

correlations,  most  of  which  are  included  In  Gough's  work.  Some  improvements 
and  changes  reflect  the  difference  of  our  approach. 

Even  with  the  reduction  of  the  problem  size  by  the  change  from  local  to 
average  functions,  one  is  faced  with  a  formidable  numerical  problem. 
Typically,  in  a  two-phase  flow  one  has  a  set  of  eleven  non-linear  partial 
differential  equations.  (Up  to  thirteen  equations  if  a  turbulence  model  is 
included ) .  In  order  to  describe  the  three-dimensional  flow  in  reasonable 
detail  one  has  to  specify  the  eleven  variables  at  a  minimum  of  about  54,000 
grid  points.  If  the  flow  is  specialized  to  axially  symmetric,  then  the 
number  of  grid  points  may  be  reduced  to  about  1,500.  Therefore,  one  should 
exploit  the  axial  symmetry  of  the  gun  whenever  possible.  The  proper 

coordinates  for  flows  with  axial  symmetry  are  cylindrical  coordinates. 

Therefore,  we  have  listed  in  Appendix  A  all  equations  in  cylindrical 
coordinates  for  flows  that  are  independent  of  the  circumferential 
coordinate. 


2.  ANALYTICAL  BASIS 


2.1  Assumptions 

In  the  next  three  Sections  (2.2,  2.3,  and  2.4)  we  shall  discuss  some 
properties  of  averaged  functions  and  develop  general  formulas  that  are 
needed  for  the  derivations  in  Section  3.  The  averages  to  be  discussed  are 
weighted  space  averages  over  a  finite  averaging  volume.  We  do  not  try  to 
establish  general  properties  of  such  averages  but  rather  concentrate  on  what 
is  needed  for  a  specific  interior  ballistics  modeling.  For  that 
application,  the  quantities  to  be  averaged  are  the  local  properties  of  a  gas 
and  of  propellant  particles  within  the  averaging  volume.  We  assume  that  no 
other  material  is  present  in  the  tube. 

The  gas  is  assumed  to  be  non— reacting  and  obeying  a  set  of  algebraic 
equations  of  state  that  permits  one  to  express  all  thermodynamic  variables 
in  terms  of  two  such  quantities.  The  particular  set  of  equations  of  state 
considered  are  the  Noble— Abel  equation  and  a  constant  ratio  of  specific  heats. 
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However,  most  of  the  results  are  independent  of  the  particular  set  of 
equations  of  state  chosen. 

We  will  assume  that  the  gas  is  viscous  and  in  a  state  without  shocks 
within  the  averaging  volume.  This  is  necessary  to  have  average  equations 
with  the  proper  differentiability  conditions.  Particular  differentiability 
conditions  of  the  local  gas  properties  will  be  enumerated  in  Section  2.2. 

If  shocks  are  present  in  the  gas  flow,  then  one  could  average  only  over 
the  shock  free  regions  and  treat  the  shocks  as  explicit  boundaries. 
However,  this  approach  has  serious  drawbacks  because  of  the  uncertainty  of 
the  corresponding  boundary  conditions  (see  Section  4.6).  Space  or  time 
averaging  is  not  the  appropriate  technique  for  the  treatment  of  interior 
ballistics  flows  with  shocks  or  other  internal  discontinuities. 

The  propellant  particles  are  assumed  to  be  incompressible  and 
elastic.  We  shall  neglect  all  effects  of  the  rotation  of  the  solid 
particles,  and  shall  assume  that  the  grains  do  not  fracture.  like  in  the 
gas,  the  local  material  properties  within  and  on  the  surface  of  each 
particle  are  assumed  to  be  differentiable  functions  of  time  and  space. 
Particulars  of  the  differentiability  conditions  will  be  enumerated  in 
Section  2.2. 


2.2  Averaging  Integrals  and  Their  Derivatives 

2.2.1  Averaging  Volume  Integrals.  We  define  the  averaging  volume  V(x) 
as  the  inside  of  a  closed  surface  S(x) .  Both  are  independent  of  time  and 
dependent  on  a  spacial  coordinate  vector  x  as  a  parameter.  For  instance,  if 
V(x)  is  a  sphere,  then  x  may  be  chosen  as  the  center  of  the  sphere.  About 
the  surface  S(x),  we  assume  that  it  has  a  well  defined  normal  almost 
everywhere.  The  shape  and  the  size  of  the  averaging  volume  are  assumed  to 
be  constant. 

The  particles  are  defined  by  corresponding  surfaces,  s  Because  the 

particles  are  moving  and  burning,  the  s^^  are  functions  of  time,  but  they 

are  independent  of  the  parametric  coordinate  vector  x.  We  assume  that  the 

particle  surfaces,  too,  have  well  defined  normals  almost  everywhere.  We 

define  as  Sp  the  union  of  all  those  particle  surfaces  spi  that  are  within 

the  averaging  volume  V,  including  its  surface  Sv.  Accordingly,  the 

intersection  S  fl  S  can  have  a  finite  area.  Most  often,  the  area  of  the 
p  v 

intersection  will  be  zero  (Figure  1). 

All  averages  will  be  defined  by  integrals  over  the  space  occupied 
either  by  gas  or  by  particles.  In  order  to  have  a  convenient  notation  for 
the  corresponding  integrals,  we  define  a  phasic  function  g  as  follows 
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Sc  ct •  2*2.1 


FINITE  AREA 
INTERSECTION 
OF  Sv  AND  SP 


ZERO  AREA 
INTERSECTION 


OF 


AND  S 


Figure  1.  Averaging  Volume 


e(t,o  = 


O  if  5  is  inside  a  particle  at  time  t  (2.1) 

1  if  5  is  outside  particles  or  on  a  particle  surface  at 
time  t. 


We  will  also  use  a  non-negative  weight  function  g  for  the  calculation  of 
averages.  Let 


VG  =  /  g(£-x)  dV(£)  »  constant  (2.2) 

V(x) 

be  the  integral  of  the  weight  function  ("the  weighted  averaging  volume"). 
Then  the  weighted  volume  fraction  occupied  by  gas  is 


ot(t,x)  = 


VG 


gas 


/  g(s-x)  dv(0 

(t  ,x) 


_L_ 

VG 


/  e(t,o  g(c-x)  dv(?) 

V(x)  (2.3) 
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Section  2.2.2 

Ihe  intrinsic  average  <j)(t,x)  of  a  function  <}>(t,x)  that  is  defined  in  the 
regions  occupied  by  gas  is  defined  by 


a(t,x)  -}>(t,x)  =  /  g(£-x)  <j>(t,£)  dV(£) 

V  Lr 

gas(t ,x) 


(2.4) 

=  ix  J  B(t,5)  g(e-x)  $(t,£)  dV(C) 

VG  V(x) 


Notice  that,  whereas  <j>(t,x)  is  defined  only  within  regions  occupied  by  gas, 
the  average  <J>(t,x)  is  defined  for  all  values  of  x  (within  limits  outlined  in 
Section  2.3). 

*  * 

A  corresponding  average  <J>(t,x)  of  a  function  <J>(t,x)  that  is  defined 
only  inside  the  particles  is  given  by 

l  l-a  (t  ,x)  ]<j)  (t  ,x)  =  /  [l-0(t,5)  ]  g  ( £  —x  )4>(t,C)dV(C)  .  (2.5) 

VG  V(x) 

Sufficient  conditions  for  the  existence  of  the  average  function  are  the^ 
piecewise  continuity  with  respect  to  x  of  the  functions  <j>(t,x)  and  <j>(t,x) 
within  their  regions  of  definition.  Obviously,  the  average  of  any  function 
of  time  only  is  the  function  itself. 

2.2.2  Time  Derivative  of  Volume  Integrals.  The  averaging  integrals 
(2.3),  (2.4),  and  (2.5)  define  functions  of  t  and  x.  In  this  section  we 

fornulate  differentiability  conditions  of  the  average  functions  with  respect 
to  time  t. 

Applying  Leibnitz  formula  (Truesdell  and  Toupin)^  to  an  averaging 

integral  (2.4)  over  V  we  obtain 

gas 


y 

C.  Truesdell  and  R.  Toupin "The  Classical  Field  Theories ,  "  in  Ens^2lop£d£a. 
o£_  PhifsicSj  S.  Flugge ,  ed. ,  Vol .  III/l,  Springer-Verlag ,  1960. 
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Sect.  2.2.2 


9 _ 

9  t 


I 


V 

gas 


(t  ,x) 


Mt.x.O  dV(0 


V 


/ 

(t,x 


9_ 
9 1 


[iKt.x.O]  dV(5) 


+  /  [♦(  t,x,5(?))(u  «n  )]  dS(c) 

Sp(t.x)  Sp  sp  (2.6) 

or 


f  6(t,0  Mt.x.o  dv(0 

9t  V(x) 


=  /  6  dV(5 )  +  /  i|»(t,x,§(?))  (u  *n  )  dS(?) 

V(x)  9t  Sp(t  ,x)  SP  SP 

where  is  the  velocity  of  a  point  of  Sp  and  ngp  is  the  outward  unit 

normal  oF  Sp  at  the  same  point.  (The  "outward"  normal  points  by  definition 
into  the  grains,  Figure  1.)  The  surface  integral  is  only  over  Sp  and  not 
over  Sv  because  the  latter  surface  is  assumed  to  be  stationary. 

The  first  integral  on  the  right-hand  side  of  Eq.  (2.6)  exists  and  is  a 
continuous  function  of  x  and  t  if  /3 1  is  a  continuous  function  of  x  and  t 
and  a  piecewise  continuous  function  of  £  .  The  surface  integral  over  Sp  in 
Eq.  (2.6)  exists  if  the  surface  velocity  is  finite.  However,  the  area  of 
the  surface  Sp  has  discontinuities  with  respect  to  x  and,  possibly,  with 
respect  to  t,  whenever  the  intersection  Sp  fl  Sv  has  a  finite  area. 
Therefore,  the  surface  integral  is  a  continuous  function  of  x  and  t  only  if 
f  =  0  on  Sv. 

Because  in  our  case 


\|;(t,x,0  =  g(?-x)  J(t,5) 


(2.7) 


we  may  formulate  the  following  sufficient  conditions  for  the  continuity  of 
the  time  derivative  of  the  averaging  integral  in  terms  of  g  and  : 


3?(t,g) 

3  t 


g(5“x) 
g(5-x)  =  0 


is  continuous  with  respect 
to  t  and  piecewise  continuous 
with  respect  to  5  in  the  domain 
of  definition  of  ^  ? 

is  continuous  in  V, 

on  the  surface  Sv. 


'  (2.8) 
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Sect.  2.2.3 

If  <j>  denotes  a  gas  phase  local  variable,  then  the  first  condition  in  Eq. 
(2.8)  applies  only  when  (£,t)  designates  a  point  in  the  gas  phase.  If 
<J>  denotes  a  solid  phase  local  variable,  then  the  first  condition  in  Eq. 
(2.8)  applies  only  when  (£,t)  designates  a  point  in  the  solid  phase. 

The  differentiation  formula  (2.6)  is  in  terms  of  g  and  <{> 


/  8g 
V 


at 


dv(°  =  a7  /  dv  -  /  gf  («  •  n  )  ds 

d  u  V  S  K  F 

P 


(2.9) 


The  corresponding  formula  for  functions  <j>  that  are  defined  within  the 
solid  grains  is 


J  (1-3)  g  ||  dV 


at 


/  (1-6)  gf  dV  +  /  gj 

v  s 

p 


(u  *  n  )  dS 
sp  sp 


(2.10) 


In  the  latter  formula,  the  surface  normal  n  again  points  into  the 
grains.  Because  now  we  are  integrating  over  the  inside  of  the  grains,  the 
sign  of  the  last  integral  in  Eq.  (2.10)  is  opposite  to  that  of  the 
corresponding  integral  in  Eq.  (2.9). 


2.2.3  Spacial  Derivatives  of  Volume  Integrals.  Applying  Leibnitz  type 
formula  to  an  averaging  integral  (2.4)  over  V^ag  one  obtains 


v  /  3(t,5)  *(t ,X,U  dv(0  =  /  6 

x  V(x)  V 


V  tpdV 

X 


J  dS 
s  -s  s 

v  p 


+ 


ibn  dS 
sp 


(2.11) 


Gauss  theorem  (Fulks 


8 


P* 


354)  applied  to  the  same  integration  volume  is 


/  8  v  MV 
v  4 


I  * 


s  -s 

v  p 


n  dS 
s 


+  J  + 


n  dS 
sp 


(2.12) 


*We  note  that  ip  and  <j>  could  he  scalars,  vectors ,  or  second  order  tensors. 
For  exajnple,  if  \p  is  a  vector ,  dots  signifying  the  divergence  of  i1  and  the 
dot  product  of  i<  should  he  used  in  Eq.  (2.11).  For  simplicity,  the  use  of 
dots  is  orr\itted  in  Section  2  wherever  if*  and  $  ewe  n°k  specified.  The 
understood  presence  or  absence  of  a  dot  should  be  clear  form  the  context. 


8W.  Fulks,  advanced  2nd  Ed.,  John  Wiley  and  Sons,  Inc.,  New  York, 

1969. 
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Sc ct«  2*2*3 

Subtracting  Eq.  (2.12)  from  Eq.  (2.11)  one  obtains 


V  /  0  i(i  dV  -  /  S  (V  +  V  )i|i  dV  -  /  i(m  dS  +  /  i|m  dS  .  (2.13) 

x  v  v  x  5  sp  sp  spfisv  sp 

Sufficient  conditions  for  the  continuity  of  the  right-hand  side  of  Eq. 
(2.13)  are 


(Vx  +  )  t(t,x,?)  is  continuous  with  respect  to 

x  and  t,  and  piecewise  continuous 
with  respect  to  £  in  the  domain 
of  definition  of 

*  0  on  Sv 


(2.14) 


In  our  application  we  want  some  of  the  average  functions  to  be 
differentiable  pwice  with  respect  to  the  spacial  variables.  By  a  formal 
differentiation  of  Eq.  (2.13)  we  obtain,  assuming  that  i|i  a  0  on  Sy, 


Next,  we  apply  the  formula  (2.13)  to  the  first  integral  on  the  right-hand 
side  of  Eq.  (2.}5)  obtaining 

v  /  e  (v  +  vr  )ip  dv  =  /  3  (v  +  vr )  (v  +  vr  )i|>  dv 

xJv  X  I  K  y  X  S  X  5 

(2.16) 


“  /  (V  +V_)i|m  dS  +  /  (V  +  V_)<|m  dS 

S  x  5  sp  s  ns  *  5  sp 

p  p1  1  V 


The  surface  integral  in  (2.15)  is 


V  f  tp  n  dS“f  7  f  n  dS  .  (2.17) 

x  c  SP  c  X  SP 


Sufficient  continuity  conditions  for  (2.16)  are 
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(V  +  V  )  (V  +  V  )f 
x  t,  x  E, 


(vx  +  =  0 


Sect.  2.2.3 

is  piecewise  continuous  with  respect  to 
£  ,  and  continuous  with  respect  to  t  j 

and  x  in  the  domain  of  definition  of  \|i  ,  I 

(2.18) 

on  •  I 


Sufficient  for  the  continuity  of  (2.17)  is  that 


V  is  continuous  with  respect  to  t  and 

x  and  piecewise  continuous  with  (2.19) 

respect  to  ^ 


Because  i|>(t,x,£)  =  g  (£ -x)£  (t  ,£  )  ,  we  may  express  the  continuity  conditions 
in  terms  of  g(£-x)  and  <|>(t,€)  as  follows. 

Sufficient  for  the  continuity  of  first  order  spacial  derivatives  of  the 
averaging  integral  is  that  (see  Eq.  (2.14)) 


?e?(t,e) 

is  piecewise  continuous  with  respect 
to  £  and  continuous  with  respect  to 
t  in  the  domain  of  the  definition  of  <j) 

g(5-x) 

is  continuous  in  V  , 

g(?-x)  =  0 

on  Sv  . 

The  integration  formula  (2.13)  in  terms  of  g  and  $  ,  if  the  conditions  (2.20) 
are  satisfied,  is 

J  egv:<t,(t,odv(0  =  v J  ei*dv  +  /  •  <2-21> 

v(x)  5  V  Sp 

* 

Ihe  corresponding  formula  to  (2.21)  for  functions  <|>  defined  within  the  solid 
grains  is 


f  [Hlg  v,^(t,0  dv(0  =V  J  [l-e]g  *  dV  -  /  g  ♦  n  ds  . 

V(x)  C  v(x)  Sp  sp  (2.22) 


Ihe  continuity  conditons  (2.18)  and  (2.19)  for  second  order  derivatives  are 
in  terms  of  g  and  <)>  as  follows 
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Sect-  2.2.4 


v£v£$(t,o 


V  g(5-x) 

X 


g(£-x)  =  0 


is  piecewise  continuous  with 
respect  to  £  and  continuous 
with  respect  to  t  in  the  domain 
of  definition  of  <j>  , 

is  piecewise  continuous  (This 
suffices  because  <j>  is 
continuously  differentiable  due 
to  the  first  condition,  Eq .  (2.20)). 

on  Sv 


(2.23) 


The  integration  formula  (2.15)  is,  if  these  conditions  are  satisfied, 


/  e^rvJ(t,5)dV(0  =  V  V  J  eg£dV  +  /  gV Jn  dS  -  /  (V  gHn  dS.  (2.24) 
V(x)  5  5  X  *V(x)  Sp  5  SP  Sp  6  SP 

In  summary,  If  the  weight  function  g  is  chosen  such  that  its  first 
derivatives  are  piecewise  continuous,  g  >  0  in  V,  and  g  =  0  on  Sv,  then  the 
averaging  integrals  are  continuously  differentiable  at  least  once  if  <j>  is 
differentiable,  and  at  least  twice  if  <j>  is  twice  differentiable  within  its 
region  of  definition. 

2*2.4  Averaging  Surface  Integrals.  Some  flow  properties  are  only 
defined  on  the  surface  of  the  propellant  grains,  e.g. ,  the  burning  rate,  the 
regression  distance,  and  the  surface  temperature.  The  corresponding 
averages  are  computed  by  surface  integrals. 

The  weighted  area  of  the  grain  surface  that  is  contained  in  the 
averaging  volume  is 


SG  =  /  g(s(t,c)-x)  dS(c)  ,  (2.25) 
Sp(t ,x) 

where  x  =  s(t,e)  defines  the  surface  and  £  represents  surface  coordinates. 
Contrary  to  the  weighted  averaging  volume  VG,  the  weighted  surface  area  SG 
is  not  a  constant  but  a  function  of  t  and  x. 


Average  surface  functions  are  defined  by 


<t>  (t  ,x) 


SG 


/  g(s(t,c)-x)4>(t,5)  dS(C) 

Sp(t ,x) 


(2.26) 
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Sect.  2. 2. A 

We  discuss  the  differentiability  of  the  surface  averages  by  considering 
a  single  grain.  Let  its  surface  s(t,?)  be  defined  in  Cartesian  coordinates 
by 


s(t  ,5  ) 


<t,5)\ 
<t,C) 
(t  ,OJ 


(2.27) 


Then  the  surface 


element  dS(?)  is  defined  by  (Courant  and  John) 


9 


dS  =  Z(t,e)d5  ,  (2.28) 

where  dc  isTthe  product  of  the  differentials  of  the  components  of  Z(t,?) 

=  (det  [(— )  (^-)])1//2,  and  3s/3c  is  the  Jacobian  matrix  of  the  function 

s(t  ,£  )  . 

The  contributions  of  the  single  grain  to  the  weighted  grain  surface 
area  is  according  to  Eq.  (2.25) 

e2 

SG.  "  //  g(s(t,c)-x)  Z(t,?)  d?  .  (2.29) 

’ 

The  time  derivative  of  SG^  is 

It  (SG1>  ■  /  (-7xS)'|f  2  dC  +  /  *  §f  «  + 

s  s 

The  integral  }.n  the  last  term  in  Eq.  (2.30)  is  to  be  taken  over  the 
intersection  C  of  the  grain  surface  s  with  the  boundary  Sv  of  the  averaging 
volume.  If  we  assume  that  g  =  0  on  Sy  then  the  integral  is  identically 
zero,  and  we  dq  not  have  to  specify  conditions  for  3C/3t. 


[  /  gZ  dC]  .  (2.30) 

Sns 

V 


9R.  Counant  and  F.  John ,  Introduction  to.  Caloidm  and  Analysis,  Vol.  II,  pp. 
459-462,  John  Wiley  and  Sons,  Inc.,  New  York,  1974. 
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Sect.  2.2.4 


Sufficient  conditions  for  the  right-hand  side  of  Eq.  (2.30)  to  be  a 
continuous  function  of  x  and  t  are 


a2s 

dzdt 

is  piecewise  continuous  with  respect  to  £ 
and  continuous  with  respect  to  t  , 

v  g 

X 

is  continuous,  with  possible  exception  of 
isolated  singular  points  , 

►  (2.31) 

g  =  o 

on  Sv 

V  g  *  0 

X 

on  Sv 

The  first  condition  in  Eq .  (2.31)  is  satisfied  if  the  grain  surface  has  a 

normal  almost  everywhere.  The  next  two  conditions  on  g(£-x)  are  essentially 
the  same  as  encountered  before  in  the  discussion  of  volume  averages.  The 

last  condition  on  g  is  new,  and  it  needs  to  be  introduced  if  Bs/3t  is  not 
equal  to  zero  and  the  intersection  sf]  $v  has  a  finite  area.  (See  the 
comment  to  Eq .  (2.6).) 

Next,  we  consider  the  spacial  derivatives  of  SG^. 
according  to  Leibnitz  type  rule 

One  obtains 

V  (SG.)  = 
x  1 

/  V  g  Z  dc  +  [  J  g  Z  dC]  ||  . 

(2.32) 

P 


The  right-hand  side  of  Eq.  (2.32)  obviously  ic  continuous  if  the  conditions 
(2.31)  are  satisfied. 

If  the  averaging  volume  contains  several  grains  then  SG  is  the  sum  of 
the  individual  SG^.  The  sum  is  continuously  differentiable  if  each  of  the 
grains  satisfies  the  first  condition  in  Eq.  (2.31) ,  and  g  satisfies  the 
other  three  conditions. 

We  now  turn  to  the  surface  average  function  $(t,x),  defined  by  Eq. 
(2.26).  We  notice  that  <|>  is  a  continuous  function  of  all  its  arguments,  if 
the  conditions)  (2.31)  are  satisifed  and  the  surface  function  4>(t,£)  is 
continuous  witl)  respect  to  time  and  piecewise  continuous  with  respect  to 
£  .  We  assume  that  £  possesses  these  properties  and  reformulate  Eq.  (2.26) 
as  follows 


♦  l  <SG.)  =  l  (f  g  ?  dS) 

i  1  is. 

Pi 


(2.33) 
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The  time  derivative  of  the  left-hand  side  of  Eq.  (2.33)  is 


Sect.  2.2.4 


\  ~  +  Ft  (SG)  +  (SG) 


(2.34) 


The  first  term  in  this  expression  is  continuous  under  our  assumptions. 
Therefore,  also  the  second  term  (and  3(|>/3t)  is  continuous,  if  the  time 
derivative  of  the  right-hand  side  of  Eq.  (2.33)  is  continuous.  The 
contribution  of  each  terra  on  the  right-hand  side  of  Eq .  (2.33)  to  the  time 

derivative  is,  vj.a  Eq.  (2.30) 


R 


ti-i  zds  +f  *!tzdc 


Pi 


Pi 


+  f  ff  +  [  J  g  ♦  Z  dC]  f£ 

Spl  »  8p^Sv 


(2.35) 


Rt^  is  a  contiguous  function  of  x  and  t  if  in  addition  to  the  condition 
(2.31)  <j)  also  satisfies  the  condition 


is  a  continuous  function  of  t  and  a 

piecewise  continuous  function  of  £  (2.36) 

on  each  Sp^. 

Because  4>  and  ($G)t,  in  Eq.  (2.34),  are  continuous  functions  if  (2.31)  and 
(2.36)  are  satisfied,  these  conditions  are  sufficient  to  insure  that  <J>(t,x) 
is  continuously  differentiable  with  respect  to  time. 

In  order  to  investigate  the  spacial  differentiability  of  <{>(t,x)  we 
differentiate  Eq.  (2.33)  with  respect  to  x.  On  the  left-hand  side  we  obtain 


3  <}) 

3  t 


\  =  4>V  (SQ)  +  ( SG)V  4> 
x  x 


(2.37) 


On  the  right-haqd  side  of  Eq.  (2.33),  each  summand  produces  the  expression 
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Sect.  2.2.5 


Rxi  -  /  (v  g)  y  z  dc  +  [  /  gf  Z  dC]V  C 


Pi 


s  .  n  s 

pi  1  V 


(2.38) 


is  continuous  if  the  conditions  (2.31)  are  satisfied.  Because 
4>V  ( SG)  is  continuous,  the  conditions  are  sufficient  for  continuous 
differentiability  of  <J>  with  respect  to  the  special  coordinate. 


Second  order  spacial  derivatives  of  surface  averaged  quantities  do  not 
enter  the  governing  equations.  Therefore,  we  do  not  formulate  existence 
condition  for  these  derivatives. 


2.2.5  Differential  Equation  for  9urface  Averages.  All  surface 
averages  satisfy  a  differential  equation  for  material  properties.  We  shall 
derive  the  equation  in  this  section. 

Let  U(t,x)  be  an  arbitrary  velocity  vector  and  let  g  satisfy  the 
conditions  (2.31).  Then  one  can  combine  Eqs.  (2.30)  and  (2.32)  obtaining 
for  the  sum  SG  pf  all  individual  SG^. 


ft  (SG)  +  yvx  (SG)  -  /(Vxg).(U  -||)  Z  dc  +/  gffdc  •  (2.39) 

S  s 

p  p 

The  integrals  on  the  right-hand  side  are  taken  over  Sp,  i.e. ,  over  all  grain 
surfaces  contained  in  the  averaging  volume. 

A  corresponding  formula  can  be  derived  for  the  product  (SG)  <j>  from  Eqs. 
(2.34),  (2.35),  (2.37),  and  (2.38)  with  the  result 


((SG)  +  U  Vx((SG)  +)  -  /  (Vxg).(U  -  $  Z  d? 

s 

p 

(2.40) 

+  /  g  ♦  ||  d£  +  /  g  ||  2  d5 

O  b 

P  P 

Next,  we  eliminate  the  derivatives  of  SG  between  Eqs.  (2.39)  and  (2.40), 
obtaining  the  differential  equation 


i!+uv<j>=—  /  gM-zdc- 

3 1  xv  SG  g  8  3  t  ^ 

P 

_ _ _ _ +-—  /  (6  -  <|>)g  — 

SG  T  f  v'fe3t 


~  /  (♦  -  ♦)  (U  -  |~)  •  (V^g)  Z  d £ 

s 

p 

(2.41) 
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The  first  integral  on  the  right-hand  6ide  of  Eq.  (2.41)  is  by  definition  the 
surface  average  of  3$/3t.  The  other  two  integrals  are  generally  assumed  to 
be  small  and  neglected  for  interior  ballistics  problems.  We  notice  that 
both  integrals  vanish  if  $  ■  on  the  propellant  surface,  i.e.,  if  the 
property  £  is  identical  for  all  grains.  If  U  is  taken  as  the  average  grain 
velocity,  the  term  U-3s/3t  may  be  small,  e.g.,  if  all  grains  have  the  same 
velocity  and  d0  not  burn,  because  3s/3t  is  equal  to  the  sum  of  the  local 
grain  velocity  and  local  surface  regression  velocity.  The  term  3Z/3t  is 
zero  if  the  grains  are  not  burning. 

If  we  neglect  the  last  two  integrals  in  Eq.  (2.41)  and  use  Eq.  (2.26) 
to  define 

♦  -w  {s#ds(5)  <2'42) 

P 

then  the  differential  equation,  Eq.  (2.41)  simplifies  to 


H +  u  v  -  <*> 


(2.43) 


where  <<j> >  is  a  model  for  <f>. 

For  the  velocity  U  one  chooses  the  average  grain  velocity,  assuming 
that  by  this  choice  one  of  the  neglected  terras  can  be  kept  small  while  not 
Introducing  another  dependent  variable. 


2.3  Regions  of  Definition  of  Average  Variables 

In  this  section  we  describe  regions  of  definition  of  the  average 
functions.  In  principle,  the  averaging  volume  V  can  be  of  any  shape  and 
size.  However,  in  order  to  preserve  an  axial  symmetry  of  the  averaged 
quantities,  the  volume  V,  the  weight  function  g,  and  the  reference  point  x 
associated  with  the  location  of  the  volume,  all  must  be  chosen  with  certain 
symmetry  properties.  Instead  of  trying  to  formulate  a  general  averaging 
volume  with  the  desired  properties,  we  give  two  examples  of  admissible 
averaging  volumes. 

The  simplest  example  of  an  averaging  volume  is  a  sphere  with  the 
reference  point  x  in  its  center  and  a  weight  function  that  depends  only  on 
the  distance  from  its  center.  Let  the  diameter  of  the  sphere  be  £. 

Another  example  is  an  orthogonal  circular  cylinder  with  the  reference 
point  at  its  center  and  with  an  axis  parallel  to  the  axis  of  the  gun  tube. 
To  be  specific,  we  assume  that  the  height  of  the  cylinder  is  2Z/3  if  i  is 
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the  diameter  of  the  cylinder.  In  this  example,  the  weight  function  depends 
on  the  radial  as  well  as  on  the  axial  coordinates  within  the  cylinder,  and 
the  volume  of  the  averaging  volume  is  the  same  as  that  of  the  spherical 
averaging  volume. 

In  both  examples,  the  quantity  £  is  equal  to  a  diameter  of  the 
averaging  volume.  In  general,  we  may  assume  a  characteristic  length 
i  associated  wiph  any  particular  averaging  volume.  The  size  of  the  volume 
and,  therefore,  the  size  of  £  ,  is  restricted  by  two  requirements.  First, 
the  averaging  volume  must  fit  inside  the  gun  barrel  and,  second,  we  want  it 
to  be  larger  than  the  largest  grain  ip  order  to  insure  that  gas  is  present 
within  every  averaging  volume.  let  D  be  the  largest  diameter  of  a  grain 
and  let  Dgun  be  the  inner  diameter  of  the  gun  tube.  Ihen  in  the  two 
examples  Z  must  satisfy  the  conditions 

^p^max  ^  ^  ^  (Dgun^min  •  (2.44) 


ftie  would  obtain  similar  restrictions  for  the  characteristic  length  of  any 
averaging  volume.  We  assume  that  D  and  D„un  are  such  that  the  inequalities 
in  Eq .  (2.44)  can  be  satisfied  by  a  margin  if  Z  is  properly  chosen. 

Hie  position  of  the  averaging  volume  (and  its  reference  point)  inside 
the  gun  tube  is  restricted.  If  a  constant  averaging  volume  intersects  a 
boundary,  then  the  sum  of  the  gas  volume  fraction  a,  as  defined  by  Eq . 
(2.3),  and  of  the  corresponding  particle  volume  fraction  is  not  equal  to 
one.  Consequently,  the  definition  of  averages  by  Eqs.  (2.2)  through  (2.5) 
cannot  be  used  if  a  non-zero  intersection  occurs,  and  the  location  of  the 
averaging  volume  is  restricted  to  positions  without  intersections  between 
the  averaging  volume  and  boundaries.  (See  also  Section  4.6)  This  means 
that  the  reference  point  x  cannot  be  moved  arbitrarily  close  to  all 
boundaries.  If  the  averaging  volume  is  a  sphere  with  the  diameter  Z ,  then  x 
is  restricted  to  locations  that  are  at  least  Z/2  away  from  the  breech,  the 
walls,  and  prqjectile  base.  In  the  second  example  (cylinder),  x  may  be 
located  at  points  that  are  at  least  Z/ 2  away  from  the  tube  walls  and  Z/3 
away  from  the  breech  and  from  the  projectile  base.  Obnsequently ,  because  of 
the  finite  sizp  of  the  averaging  volume,  none  of  the  averaged  quantities  are 
defined  in  the  boundary  regions.  If  the  grain  diameter  D  is  large,  then 
the  regions  where  the  averaged  quantities  are  not  defined  can  be  a 
significant  part  of  the  interior  of  the  gun  tube. 

In  the  remaining  regions,  the  porosity  a  and  all  averages  pertaining  to 
gas  properties  are  everywhere  defined  by  Eqs.  (2.3)  and  (2.4),  respectively. 

Average  properties  of  propellant  grains  are  defined  by  Eq .  (2.5).  The 
definition  provides  a  value  for  the  average  function  only  if  a  <  1,  i.e. ,  if 
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there  are  grains  within  the  averaging  volume.  The  limitation  also  holds  for 
surface  averaged  quantities,  defined  by  Eq.  (2.26).  The  surface  averaged 
quantities  are  grain  properties  and  they  are  defined  only  if  there  are 
grains  within  the  averaging  volume. 

Another  average  dependent  variable  which  is  introduced  in  Section  4.2 
is  the  weighted  number  m  of  grains  in  the  averaging  volume  that  is  defined 
by 


m(t,x)  =  VG  (l-a)/v  (d)  ,  (2.45) 

P 

A  * 

where  d  is  the  average  regression  distance  of  the  grains  and  v  (d)  is  the 
corresponding  grain  volume,  given  by  a  correlation  function.  According  to 
the  definition,  m  is  indeterminate  in  regions  without  grains,  because  d  is 
not  defined  in  those  regions.  We  notice,  however  that  m  >  0  and  V^m  -*  0  as 
x  moves  to  a  position  where  the  averaging  volume  contains  no  grains. 
Therefore,  we  may  define  a  continuation  m  =  0  in  regions  without  grains. 
With  this  extension,  m  is  defined  in  all  those  regions  where  gas  properties 
are  defined,  i.e.  ,  everywhere,  except  in  boundary  regions. 


2.4  Averaging  Weight  Function 

The  averaging  weight  function  g(y)  is  defined  inside  the  averaging 
volume  V  and  on  its  boundary  Sy.  It  has  the  following  properties  (see 
Sections  2.2.2,  2.2.3,  and  2.2.4) 


g  >  0 

in  V 

g  =  0 

on  Sv  , 

(2.46) 

Vg 

continuous  in  V  with  possible  exception  of  isolated 
singular  points  , 

Vg  =  0 

on  S 

V 

Next,  we 

give  examples  of  functions  g(y)  that  satisfy  these  conditions 

the  two 

examples  of  averaging  volumes  mentioned  in  the 

previous 

section.  Let  y  =£-x,  i.e.,  let  the  point  of  origin  of  the  coordinate 

vector  y  be  at  the  center  of  the  averaging  volume.  (In  both  our  examples 
the  center  coincides  with  the  reference  point  x.) 
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If  V  is  a  sphere  with  the  diameter  £,  then  we  may  define  the  weight 
function  by 


g(y)  =  (2+n)  (3+n)  (4-fn)  ^  ^ 


-  11 


1+n 


for  -  j  <  y  <  Z/2 


(2.47) 


with  an  arbitrary  n  >  0.  The  weighted  averaging  volume  VG  is  for  this  g(y) 

Z/2  ?  4  n  o 

VG  =  /  g  dV  =  4it  /  g(y)y  dy  =  -^  tt  (y)  .  (2.48) 

V  o 

As  a  second  example  we  chose  a  cylinder  with  the  diameter  Z  and  height 
2Z/3*  Let  r  and  z  be  the  radial  and  axial  coordinates  within  the  cylinder, 
with  the  point  of  origin  at  the  center  of  the  cylinder  Then  we  may  define 
with  arbitrary  positive  m  and  n 


g(r,z)  ={(2+m)  (2+n)  (3+n)  (  1  -  1^)  1_hn(  1  -  i^)  1+n  .  (2.49) 

The  weighted  averaging  volume  VG  is  for  this  choice  of  g 
Z/3  Z /2 

VG  =  /  g  dV  =  4ir  J  /  g(r,z)r  drdz  =  ^  Tt(y)  3  ,  (2.50) 

V  o  o 

i.e.,  equal  to  the  volume  |v|  of  the  cylinder  itself. 

In  both  examples,  we  have  weight  functions  with  a  maximum  at  the  center 
of  the  averaging  volume.  The  functions  are  continuous  but  their  gradients 
possess  discontinuities.  The  weight  function  for  the  spherical  averaging 
volume  has  a  discontinuous  point  at  the  center  of  the  sphere.  The  second 
weight  function  has  a  singular  gradient  along  the  line  r  =  0  and  on  the 
plane  z  =  0.  Therefore,  if  the  flow  includes  phenomena  that  require  surface 
averaging  one  should  use  a  different  weight  function  for  the  cylindrical 
averaging  volume.  (For  volume  averaging,  piecewise  continuity  of  Vg  is 
sufficient . ) 

The  following  two  weight  functions  have  no  discontinuities.  They  are 
chosen  such  that  the  weighted  averaging  volume  is  the  same  as  before,  i.e., 
equal  to  the  volume  of  a  sphere  with  diameter  Z . 
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A  weight  function  example  for  a  sphere  is 


(2.51) 


A  weight  function  for  the  cylindrical  averaging  volume  is 


g(r,z)  =  —z — 
TT  -4 


[cos  (jyj)  +  lj  [cos  +  1] 


(2.52) 


Numerous  other  examples  can  be  constructed,  e.g. ,  based  on  the 
functions 


g(r)  -  (1  -  (^72)2m)1+n 


(2.53) 


(2.54) 


and  corresponding  for  the  dependence  on  z.  Particularly,  functions  of  the 
type  (2.53)  with  large  integer  m  and  small  positive  n  have  properties  that 
are  desirable  according  to  Section  4.2.1. 


3.  CONSERVATION  EQUATIONS 


The  mathematical  description  of  a  two-phase  flow  field  is  composed  of 
two  sets  of  local  conservation  equations  (one  for  each  phase)  ,  a  set  of 
local  constitutive  relations  for  each  phase,  and  interfacial  or  jump 
conditions  which  relate  locally  the  two  phases  only  on  the  interfaces.  As 
in  other  two— phase  models  of  interior  ballistics,  all  chemical  reactions  are 
excluded.  Burning  of  the  grains  is  represented  by  a  transfer  of  mass, 
momentum,  and  energy  from  the  solid  phase  to  the  gas  phase.  Furthermore, 
the  effects  of  body  forces  on  both  phases  are  assumed  to  be  negligible.  By 
averaging  the  local  conservation  equations  according  to  the  definitions  and 
formulas  determined  in  Section  2,  and  by  using  the  local  interfacial 
conditions,  we  derive  the  coupled  set  of  average  two-phase  equations.  The 
details  of  this  procedure  are  given  in  this  section.  The  average  equations 
in  vector  form  are  derived  in  three  spatial  dimensions  and  time.  The 
governing  equations  for  axially  symmetric  flow  in  cylindrical  coordinates 
are  listed  in  Appendix  A. 


25 


3. 1  local  Equations 


Sect.  3.  1.  1 


3.1.1  local  Conservation  Equations.  The  flow  field  is  assumed  to  be 
composed  of  two  disjoint  phases:  gas  and  solid  grains.  The  gas  is  assumed 
to  be  compressible,  viscous  and  heat  conducting.  The  local  conservation 
equations  for  the  gas  are  the  Navier- Stokes  equations  (Tsien,  pp,  3—16) 


3“  +  V*  (pu)  =  0  ,  (3.  1) 


8  (pue) 
8t 


V  vv 

+  V*  (puu) 


-  Vp  +  7»n 


(3.2) 


8  (pe) 

at 


+  V*  (pue)  -  - 


p  V*u  +  }  -  V'Q 


(3.3) 


where  p,  e,  and  u  are  the  density,  specific  internal  energy,  and  the 
velocity  vector,  respectively.  The  constitutive  laws  for  the  viscous  stress 
tensor  II,  the  heat  dissipation  function  and  the  heat  conduction  vector 

Q  are 


ff  =  2y E  +  (X  -|y)  V*u  I  ,  (3.4) 

N/  X/  X/  X#  ^  X/  9  ^  O 

*j  =  2vi  E:E  +  (X  -  |  y)  (V-u)Z  ,  (3.5) 


Q  =  -  k  VT  ,  (3.6) 


where 


E  =  0.5  [Vu  +  (Vu)TJ 


(3.7) 


V  V  V 

and  y,  X,  <  are  the  shear  viscosity  coefficient,  the  bulk  viscosity 
coefficient  and  the  heat  conduction  coefficient,  respectively,  that  may 
depend  on  the  local  temperature  T.  The  local  pressure  and  temperature  are 
given  by  equations  of  state  of  the  form  p  =  p(p ,e)  and  T  =  T(p,e)  . 


Tsien,  " The  Equations  of  Gas  Dynamics ,  "  in  Fundamentals  of  Gas 
Dynamics,  "  H*  W.  Emmons,  ed»,  Princeton  University  Press,  Princeton,  NJ t 
UW. 
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Each  solid  grain  is  assumed  to  be  incompressible  (the  density  of  a 
* 

grain  p  =  constant)  but  deformable.  The  local  conservation  equations  for 
the  solid  phase  can  be  expressed  in  a  form  similar  to  those  of  Eqs.  (3.1) 
and  (3.2)  (Prager) 

(p)  +  V.  (p  u)  =  0  ,  (3.8) 


9  *  *  *  *  *  * 

(p  u)  +  v*  (p  u  u)  =  v*n 


(3.9) 


where  u  is  the  local  velocity  vector  of  the  grain.  For  our  purposes,  the 

: k 

solid  phase  stress  tensor  H  represents  the  total  stress  within  the  solid 

* 

grain.  A  constitutive  law  for  II  could  be  based  on  Hooke’s  law.  Although 
the  local  angular  momentum  of  the  grains  could  be  significant,  it  is  assumed 
that  the  average  effect  of  the  angular  momentum  is  small  and  can  be 
neglected.  Cbnsequently ,  the  local  conservation  equation  for  the  angular 
momentum  of  a  grain  is  omitted. 

3.1.2  local  Interfacial  Conditions.  The  interfacial  conditions  relate 
the  two  disjoint  phases.  The  interface  between  the  gas  and  solid  is 
considered  a  singular  surface  across  which  mass,  momentum  and  energy  is 
transferred.  The  conditions  tha^  are  valid  on  the  interface  can  be 
expressed  as  (Truesdell  and  Toupin  )  : 


*  * 


n*p  (u  -  u  )  =  n*p(u  -  u  )  , 

sp  sp 


n*  p  (u  -  u  )u  +  np-n*n=n*p(u-u  )  u  -  n»  II 


sp 


sp 


(3.10) 

(3.11) 


^  i  ^  ^  ^  ^  ^  ^  ^ 

n* p  (u  -  u  )  te  +  y  u*  u]  +  P  n»  u  +  n»  Q  -  n* n  •  u 


A  A  V  *  *  *  £ 

=  n  •  p  (u  -  u  )  [e+-^-u*u]+n*Q-n*n*u 
sp  Z 


(3.12) 


11 W.  Prayer,  Introduction  to  Mechanics  of  Continua,  Ginn  and  Company ,  New 
York,  1961. 
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X 

where  u  is  the  local  interface  velocity,  n  is  a  unit  normal,  and  Q  is  the 
sp 

local  heat  conduction  vector  within  the  grain. 

The  local  interface  velocity  u  is  defined  in  terms  of  the  local 

sp 

regression  rate  d  of  the  grain  surface 

v  V 

u  ft, 5(0)  -  S(t,5(?))  +n  d(t,5(0)  ,  (3.13) 

s  p  s  p 

where  5  is  the  surface  coordinate  vector,  d  >  0  and  n  is  the  unit  normal 

s  p 

to  the  grain,  outward  with  respect  to  the  gas. 


3. 2  Averaging  of  the  Local  Conservation  Equations 

3.2.1  Derivations  of  the  Average  Gas  Continuity  Equation  and  Porosity 
Equation.  To  derive  the  average  gas  phase  continuity  equation,  we  multiply 
Eq.  (3.1)  by  3  (t  ,5 )g (5 -x) ,  integrate  over  the  averaging  volume  V(x)  and 
obtain 

/  g(t,5)g(5-x)  3.p(l;;g)  dv(5) 

V(x)  9t 

(3.14) 


+  /  3  (t  ,5)g(5“x)  V  •  [p  (t  ,5  )u(t  ,5)  J  dV(5)  =0 

V(x)  * 

Using  formulas  (2.9)  and  (2.21)  with  respect  to  the  first  and  second 
integrals  of  (3.14),  respectively,  we  have 


3  ^  ^  ^ 

r—  /  3  gp  dV  +  V  •  /  3gpu  dV  +  /  g  p  (u  -  u  )*n  dS=0  .  ((3.15) 

0  1  y  X  ^  g'  SP  SP 

P 

By  the  definition  of  a  volume  averaged  quantity  (2.4)  and  the  interfacial 
mass  flux  condition  (3.10),  Eq.  (3.15)  can  be  written  as 

|-£-  (a  (t  ,x)p  (t  ,x))  +  V*  (a(t,x)[pu]  (t,x)) 


+ 


* 

P 

VG 


.  * 
/  g(u 
S 


P 


u  )•  n  dS  =  0 
sp  sp 


(3.16) 
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=  constant  and  VG  =  constant.  In  Eq.  (3.16)  p  is  the  average 

is  the  average  of  the  gas  momentum  density 


because 

gas  density  and  the  quantity 
u.  We  define  the  average  gas  velocity  vector  u  as  the  ratio 


pu 


u(t ,x) 


pu 

(t  ,x) 

p< 

.t  ,x) 

(3.17) 


Using  this  definition  of  u,  the  local  regression  rate,  defined  by  Eq . 
(3.13),  and  the  definition  of  the  average  surface  function  (2.26),  we  can 
rewrite  the  average  gas  continuity  Eq.  (3.16)  as 

[a(t  ,x)p  (t  ,x)  ]  +  V-  [a  (t  ,x)p  (t  ,x)u(t,x)  ]  =  p  SG%-~  d(t,x)  .  (3.18) 

The  derivation  of  the  average  solid  phase  continuity  equation  proceeds 
in  a  similar  fashion  to  that  of  the  average  gas  continuity  equation. 
Multiplying  Eq.  (3.8)  by  (l-g)g,  integrating  over  V(x) ,  invoking  formulas 
(2.10)  and  (2.22),  and  using  the  definition  (2.5)  of  an  average  solid  grain 

property,  and  (3.13)  of  the  local  regression  rate,  we  have 

(VG(l-ct)p)  +  V*(VG(l-a)  |  pu|)  -  p  /  g  d  dS  =  0  .  (3.19) 

S 

P 

* 

Using  the  surface  average  definition  (2.26)  and  the  fact  that  p  is  a 
constant,  Eq.  (3.19)  can  be  written  as 

(l-a)  +  V.[(l-a)5]  =  ||  d  .  (3.20) 

Hence,  for  incompressible  solid  grains,  the  average  continuity  equation  for 
the  solid  phase,  Eq.  (3.20),  is  the  governing  equation  for  the  porosity  a. 

We  notice  that,  if  the  density  is  constant  or  depends  only  o£i  time, 

then  the  average  velocity  is  given  directly  by  Eq.  (2.5),  e.g.,  u.  The 

different  definition  of  the  average  gas  velocity  via  the  average  momentum 
density  by  Eq .  (3.17)  is  advantageous  when  the  density  depends  on  the 
spatial  coordinate. 

The  average  gas  continuity  equation,  Eq.  (3.  18) ,  is  coupled  to  the 

A  # 

solid  phase  by  the  source  term  p(SG/VG)d.  As  expected,  the  amount  of  mass 
added  to  the  gas  phase  is  exactly  the  amount  liberated  from  the  solid 

phase.  If  the  grains  are  not  regressing  (not  burning)  ,  then  the  average 

. 

regression  rate  d  and  the  source  term  are  zero.  The  surface  average 
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SG  and  the  surface  average  regression  rate  d  are  two  new  unknowns.  To 
restrict  the  number  of  unknowns,  d  is  replaced  by  a  correlation  (denoted  by 
<d>)  which  is  obtained  from  experiments  (see  Section  4.7.7).  To  understand 
the  error  involved  in  such  a  substitution,  we  rewrite  Eq.  (3.  18)  as 


l^lapl  +  V*  [apu]  =  p  <d(t,x)> 


(3.21) 


+  S  f§  [-|q  /  d(t,5(c))  dS(c)  -  <d(t,x»]  . 

s 

p 

The  bracketed  term  on  the  right-hand  side  of  Eq.  (3.21)  is  the  error  term 
and  is  equal  to 

—  I  [  d  J  g(5(?)-x)  dS(Oj  -  <d(t,x)>  (3.22a) 

SG  i  1  s  . 

Pi 

12 

by  the  mean  value  theorem  for  multiple  integrals  (Apostol)  and  where  5^  is 
some  point  on  From  expression  (3.22a),  the  following  inequality  can  be 

derived  : 


/  g  d  dS(c)  -  <d>  I  <  max|d(  t,C(C .))  -  <d(t,x)>j  .  (3.22b) 

i SG  s  i  1 

P 

Thus,  a  sufficient  condition  for  the  error  to  be  small  is  that  the 

*  v 

difference  between  the  local  regression  rate  d  over  each  surface  and  the 
value  of  the  correlation  <d>  at  point  x  is  small.  A  common  expression  for 
<d>  is  given  by  Eq.  (4.100).  If  the  error  given  by  Eq.  (3.22a)  is  not 
small,  another  correlation  for  <d>  must  be  used.  In  practice,  the  error  is 
assumed  small  and  Eqs.  (3.18)  and  (3.20)  are  written  with  d  replaced  by 
<d>.  Furthermore,  an  additional  formal  error  could  be  introduced  by  the 
modeling  of  SG.  However,  this  is  avoided  by  the  definition  of  m  in  terms  of 
SG  (see  Section  4.7.8). 

3.2.2  Derivations  of  the  Average  Gas  and  Solid  Momentum  Equations. 

The  average  gas  momentum  equation  is  derived  by  multiplying  the  local 
momentum  equation,  Eq.  (3.2),  by  the  function  gg,  by  integrating  over  the 
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averaging  volume  V(x)  and  by  applying  formulas  (2.9)  and  (2.21).  The 
results  of  these  operations  can  be  written  as 


9_ 

3t 


JggpudV-Jgpun  *u  dS  +  V  •  /  3g  puu  dV 

v  s  sp  sp  X  V 

P 


+  J  g  n  *puu  dS  =  -  V  J  0g  p  dV  +  V  •  J  gg  n  dV 

s  sp  x  V  x  V 

p 


g  ("spp 


n  »ff)  dS 
sp 


(3.23) 


We  use  the  definition  of  an  average  gas  property  (2.4)  and  the  definition  of 
u  (3.17)  in  Eq.  (3.23)  to  obtain 


9 _ 

9 1 


[a  (t  ,x)p  (t  ,x)u(t  ,x)  ]  +  V*  [a  (t  ,x) 


p  uu 


(t ,x)  ] 


=  -  V  {^—  /  egp  dV}  +  v*  {rj—  J  egff  dv}  (3.24) 

V  V 

-  -77;  /  g{n  p  -  n  *n+n  •  p  [  u  -  u  luldS 

VG  ;  L  spr  sp  sp  L  spJ  J 

S 

P 

The  term  I  puu  I  (t  «x)  represents  the  average  of  the  tensor  puu.  Because 
the  average  quantities  p  and  u  are  already  defined,  we  can  denote  the 
fluctuations  of  the  values  of  the  local  variables  from  the  value  of  the 
average  variables  as 

p'(t,5,x)  =  p  (t  ,5)  “  P  (t  ,x)  , 


and 


(3.25) 


U*(t,5,x)  =  u(t,£)  -  u(t,x) 


If  we  substitute  fornulas  (3.25)  into  the  integral  representation  of  q  fpuu]  , 
we  obtain 


1 

VG 


.  W  V  |  .  N/  V 

J  3g  puu  dV  =  apuu  +  7777  J  3g  pu'u'  dV 
V  V°  V 


(3.26) 
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The  difference  between  the  first  terra  on  the  right-hand  side  of  Eq ,  (3.26) 
and  the  left-hand  side,  involves  a  volume  average  of  the  product 
of  velocity  fluctuations.  We  define  this  difference  as  the  turbulent  stress 
tensor  of  the  flow.  Thus,  turbulence  in  this  report  is  defined  as  volume 
averaged  fluctuations.  The  turbulent  stress  tensor  11^,  models  the  quantity 

-  —  -tr  /  3 g  pu'u'  dV  =  puu  -  |p  uu|  .  (3.27) 

a  V(j  V 

We  shall  not  discuss  particular  turbulence  models  in  this  report.  A  model 

is  proposed  in  GLbeling  et  al.^  Iking  the  integral  representation  of  ^[jTuul 

and  applying  the  mean  value  theorem  for  multiple-integrals  (Apostol),  we 

can  rewrite  the  right-hand  side  of  Eq .  (3.27)  when  is  a  connected  set 

gas 

as 


[u(t  ,x)u(t  ,x)  -  u(t  ,5  )u(t  ,5  )  ]  p(t,x) 


(3.28) 


where  £  lies  in  V___  and  is  different  for  each  component  of  the  tensor  uu. 
gas 

From  Eq.  (3.28),  a  good  model  of  the  turbulent  stress  tensor  for 
compressible  flows  is  one  which  models  the  significant  differences  between 
the  tensors  uu  and  uu.  With  respect  to  the  errors  generated  by  such  a  model 
II  T  in  Eq .  (3.24),  we  want  the  errors  in  the  vector 


V •  { all  ^  -  [apuu  -  a[  puu] } 


(3.29) 


to  be  minimized  by  the  model. 

Substituting  Eq .  (3.27)  into  Eq .  (3.24),  using  Eq .  (3.13),  and 

algebraically  manipulating  the  result,  we  have 
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*  [apu]  +  V*  [apuu  J  =  -  aVp  +  V«  (all )  +  V*  (aIIT) 
o  t  i 

+ 


+ 


+ 


+ 


where  p  and  II  are  the  conctitutive  models  for  the  average  pressure  defined 

by  [JyBgpdV] / [a* VGJ  and  the  average  viscous  stress  tensor  defined  by 

[/yftgTIdV] /(a« VGJ  ,  respectively.  In  general,  it  is  simpler  to  model  the 
average  pressure  and  viscous  stress  tensor  than  to  actually  integrate  the 
local  constitutive  laws.  Each  term  in  Eq.  (3.30)  which  is  enclosed  by 
braces  is  an  error  term.  We  now  shall  discuss  each  error. 

*  • 

The  errors  in  the  models  p,  II,  H^,  and  those  introduced  by  u<d>  are 

represented  by  the  last  four  terms  on  the  right-hand  side  of  Eq.  (3.30).  If 

V__0  is  connected,  the  errors  in  the  last  two  terms  can  be  written  as 
gas 

Vx  /  Bgp  dV  -  ap]  =v{a(t,x)  [p(t,£(x))  -p(t,x)]}  (3.31) 

and 

i  a 

V  •  [t77t  /  ggn  dV  -  an]  =  V«{a(t,x)[ff(t,£(x))  -n(t,x)]}  ,  (3.32) 

x  VLr  v 

A 

where  £  (x)  are  the  mean  value  points  in  V  _(t,x)  which,  in  general,  are 

different  for  p  and  for  each  components  of  the  tensor  II.  The  models  p  and 

II  as  well  as  the  errors  (3.31)  and  (3.32)  are  discussed  in  Sections  4.7.1 


*  CY2  -k  .  r  1  r  .  •> 

P  W  u<d>  ‘  tvc  I  8("»nP  -  "s„-")dS  +  P 


VG 


sp  sp 


fl  f  r  v  v  v  v  k  k  ^  k  __ 

1 7777  J  g  n  •  p  (u-u  )u  -  n  *p(u-u  )u  dSl 
1  VG  J  OL  sp  sp  sp  sp  J  J 

J 
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r  1  r  *  •  *  SG  *  v  i 

{ VG  /  p  d  u  dS  “  ^  u<d>l 


VG 


(3.30) 


/  gg(puunpuu)  dV  -  anT]} 


/  3gp  dV  -  ap]} 


V-[i^  /  ggn  dV  -  an]} 
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and  4.7.2,  respectively.  For  the  gas  momentum  equation,  the  best 
approximations  for  p  and  II  are  the  ones  which  minimize  both  the  differences 
in  their  values  and  their  derivatives  . 

The  error  in  the  turbulence  model  was  discussed  previously  in  this 
section . 

The  second  braced  terra  in  Eq  .  (3.30)  can  be  written  as 
p  W  h  4s  /  g  t;  (C ))  ds(c)  -  u(t,x)  <d(t,x)>}  .(3.33) 

i  3  spi 

•w  v< 

* 

If  both,  u  and  d  which  are  defined  on  the  grain  surface  are  functions  of  t 
only,  then  expression  (3.33)  is  zero  and  no  error  exists.  When  this  is  not 
the  case,  one  can  bound  (3.33)  using  the  mean  value  theorem  for  multiple 
integrals  by 

i *  sn i  i  * ,  n  •  *  •  , 

|p  — |  max  |u(t,C(CjL))  d(t,x)  -  u(t,x)  <d(t,x)>|  ,  (3.34) 

where  is  different  on  each  surface  sp^ .  Expression  (3.34)  can  be  bounded 
by 

l  ^  QfJ  lr*  |  /»  \  A  |  ll*  •  |  . 

|p  —  I  {  d(t,x)max|u(  t ,5 (C i))  -  u(t,x)|  +  | u( t ,x) | | <d( t ,x)>  -  d(t,x)|}  . 

(3.35) 


Thus,  the  error  in  replacing  -gg-  J  s^gu  d  dS  with  u<d>  consists  of  two 
parts.  One  error  involves  the  approximation  of  d  with  <d>  and  is  discussed 
in  Section  3  .2  .1  .  The  other  term  is  small  if  the  values  of  the  local 
particle  velocity  at  the  grain  surfaces  are  near  that  of  the  average 

particle  velocity  at  x;  that  is  if  the  fluctuations  are  small.  If  Jboth 

*  . 

terras  are  not  small,  then  a  correlation  of  the  fluctuations  u  dT  must  be 
modeled  and  included  in  Eq .  (3.30). 

The  terra 


}_ 

VG 


/  g  [» 


sp 
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n  • 

sp 


V  V 


*  *  V, 

p  (u-u 


sp 


)u]  dS 


(3  .36) 
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can  be  rewritten  using  the  mass  flux  jump  Eq.  (3.10)  and  regression  rate 
definition  (3.13)  as 

■  p  ^  /  g  (u  -  u)  d  dS  ,  (3.37) 

VG  S 

P 

or  using  the  momentum  flux  jump  Eq .  (3.11)  as 


h  /  g  [(nsp#fi  '  nsp  ”  nsp*^  dS  •  (3-38) 

S 

P 

On  the  interface  between  the  gas  and  the  particles,  we  assume  either  that 
the  normal  stresses  are  equal  (the  integrand  in  Eq.  (3.38)  is  zero),  or 
equivalently,  that  the  gas  and  particle  velocities  are  equal  (the  difference 

in  the  integrand  in  Eq .  (3.37)  is  zero).  In  the  special  case  of  no  burning 

• 

d  =0,  the  error  is  zero.  When  the  above  assumption  is  not  true,  the 
expression  (3.36)  must  be  modeled  by  a  correlation. 

v 

From  Eq.  (2.21)  with  <f>  =  1 ,  we  have  the  relationship 


Va 


1 

VG 


dS 


Using  the  formula  (3.39),  we  have  the  equality 


(3.39) 


h  /  glnspp _  nSp*]  ds  +  pVa 

S 

P  (3.40) 

=  h  {  gtnsP(p " p)  "  Vn3  dS  • 

b 

P 


We  define  the  surface  integral  on  the  right-hand  side  of  Eq .  (3.40)  as  the 
drag  force.  Ihe  drag  force  is  modeled  by  the  correlation  D  which  is 
discussed  in  Section  4.7.5.  The  error  incurred  by  this  approximation  is 


{  /  g[nsp(p  -  p)  -  nsp*n]  ds  -  D(t,x)} 

s 

p 


(3.41) 
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This  definition  is  consistent  with  Ishii's^  development  but  is 
different  from  Gibeling  et  al.  and  Gough' s*  which  is  defined  in  terms  of 
the  surface  integral  of  the  weighted  fluctuation  of  the  normal  total  gas 
stress  tensor;  ng  ^  (IHH )  -  ngp(p-p).  For  the  special  case  when  the  average 
viscous  stress  tensor  is  zero  (the  inviscid  two-phase  model),  our  definition 

C  1 

and  those  of  Gibeling  et  al.  and  Gough1  agree.  We  recognize  the  fact  that 
Eq.  (3.41)  is  a  formal  definition  which  may  not  correspond  to  an 
experimentally  determined  drag  force.  In  such  a  case,  other  effects 
included  in  the  experimental  drag  force  would  have  to  be  subtracted  to 
obtain  the  correlation  corresponding  to  D. 

The  derivation  of  the  average  solid  phase  momentum  equation  parallels 
that  for  the  average  gas  momentum  equation.  We  multiply  Eq .  (3.9)  by 
(1-0 )g,  integrate  over  the  averaging  volume  V,  use  formulas  (2.10)  and 
(2.22)  and  the  definition  of  the  average  of  a  solid  grain  property  (2.5)  to 
obtain 


l~i  [(1_a 


“ST 

pu 


(t,x)]  +  V*  [  ( 1  -a)) 


j  *** 

1  puu 


(t ,x)] 


=  v  •  I  (l-6)gn(t,5)  dv}  (3.42) 

A  v  u  y 

V  ^  V  x# 

+  J  g  nsp.p*(u-usp)S  dS  -  /  g  nsp.n  dS 
P  P 


Because  p  is  a  constant, 


pu 


k  k 

(t ,x)  =  p  u(t ,x)  and 


puu 


(t,x)  =  p 


TF* 

uu 


(t  ,x) 


By  adding  and  subtracting  V[(l-ct)p],  by  using  Eq.  (3.39),  and  replacing 
* 

nSp*n  on  the  surface  with  its  equivalent  via  the  momentum  flux  interfacial 
jump  condition  (3.11),  we  can  rewrite  Eq .  (3.42)  as 


15 

M.  Ishii,  Thermo-Fluid  Dynamic  Theory  of  Two-Phase  Flow,  Eyrolles,  France 
1975. 
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|^-[(l-a)pu]  +  V*  [  (l-a)p_uu  ]  =-  (l-a)Vp 


+  Vx’^vg  I  (i-e)gn  dV  +  (l^t)  pi} 


+  VG  /  g  nsp*p*(^sp^  dS 

J 

P 

+  /  gfn  p-n  p-n  *ff]  dS 

VG  J  6L  spF  sp  sp  J 


VI  V  V  V  *  *  %#  *1 

+  VTTT  J  g  n  *p  (u-u  )u  -  u  •pCu-u  )u|  dS  , 
VG  Js  &L  sp  sp  sp  sp  J 
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where  I  is  the  identity  tensor.  Eq .  (3.43)  can  be  rewritten  as 


[  (l-a)pu]  +  V*  [  (l-tx)puu]  =  -  (l-ot)Vp  +  V*  [(l-a)fi] 


***, 


r  *  ,  *  CQ  *  •  ] 

+  V»  [  (l-a)nT]  -P  <d>  +  ^  D 
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VG 
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** 
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+  {v*[(ic  l  g[nsP(p_p)  _  ds  ~h 


where  II  is  the  constitutive  model  for  the  average  stress  tensor 
solid  phase  and  represents 


I^tIg  /  «-«>*"  dv  +  P1  • 


(3.43) 


(3.44) 


for  the 


(3.45) 
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* 

and  II ^  is  the  constitutive  model^  for  the  average  solid  phase  turbulent 
stress  tensor.  In  analogy  to  n T ,  II  models  the  tensor  (see  Eq.  (3.27)) 


(t,x)  *  "T^TW  /  u'  *'  dV  .  (3.46) 

* 

We  recall  that  by  definition  II  denotes  the  total  stress  tensor  for  the  solid 

grain.  The  quantity  defined  by  Eq.  (3.45)  is  the  difference  between  the 

* 

average  total  stress  in  the  solid  phase  (the  integral  of  (l“g)gII/VG  over  the 
averaging  volume)  and  the  stress  caused  by  the  average  gas  pressure  (-pi). 
The  resulting  stress  is  the  stress  caused  by  the  grains  themselves,  for 
example,  by  the  compactif ication  of  the  propellant  bed.  Consequently,  we 
call  the  expression  (3.45),  the  average  intergranular  stress,  and  II  the 
average  intergranular  stress  model.  As  with  the  average  pressure,  viscous 
stress  tensor,  and  turbulent  gas  stress  tensor,  it  is  simpler  to  separately 
model  the  intergranular  stress,  the  solid  phase  turbulent  stress,  and  the 
drag.  The  errors  incurred  by  these  models  are  represented  by  the  last  three 
braced  terms  in  Eq.  (3.44). 

The  remaining  error  terms  in  Eq.  (3.44)  (those  enclosed  by  braces)  are 
the  surface  integral  involving  the  velocity  or  stress  jump,  and  the  surface 
integral  representing  the  source  term.  These  terms  are  discussed  in  the 
derivation  of  the  average  gas  momentum  equation  (see  the  analyses  beginning 
near  Eqs.  (3.36)  and  (3.33),  respectively). 

3.2.3  Peri vat ion  of  the  Average  Gas  Internal  Energy  Equation.  The 

average  internal  energy  is  needed  to  compute  certain  quantities,  e.g.,  the 
pressure  and  temperature  via  the  equations  of  state  for  the  average 
quantities.  The  average  internal  energy  can  be  obtained  in  either  of  two 
ways.  First,  by  adding  the  local  internal  energy  equation  to  the  equation 
for  the  local  kinetic  energy,  an  equation  for  the  local  total  energy  can  be 
written.  Following  a  similar  procedure  to  those  given  in  Sections  3.2.1  and 
3.2.2,  we  then  can  derive  an  average  total  energy  equation.  Finally,  the 
average  internal  energy  value  is  obtained  as  the  difference  between  the 
average  total  energy,  and  the  average  kinetic  energy  determined  from  the 
average  velocities.  The  second  way  is  to  average  the  local  internal  energy 
equation,  Eq.  (3.3),  directly.  The  former  procedure  is  the  most  common. 
However,  we  use  the  latter  approach  because  several  terms  which  must  be 
assumed  small  or  modeled  by  additional  correlations  can  be  avoided,  and  the 
terms  which  must  be  modeled,  have  simpler  physical  interpretations,  and 
therefore,  are  easier  to  model.  An  example  of  a  term  that  can  be  eliminated 
by  the  second  method  but  is  present  in  the  first  is 


jc  "k 

u(t,x)u(t,x)  - 


** 

Ult 
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/  S(t,5)  g(f;-x)  [p(t,5)u(t,5)*u(t,5)  -  p (t,x)u(t,x)*u(t,x)]  dV 
V 

(3.47) 


=  /  Sg  p  u’-u'  dV 
V 

The  non-negative  integral  (3.47)  is  the  average  difference  between  the  local 
kinetic  energy  and  the  dot  product  of  the  average  velocity  times  the 
density.  An  example  of  a  term  that  can  be  modeled  more  easily  in  the 
average  internal  energy  equation  is  the  dissipation  term.  In  the  average 
internal  energy  equation,  the  term  $  represents  the  average  conversion  of 
viscous  work  by  the  fluid  into  heat  only,  whereas,  in  the  average  total 
energy  equation,  the  corresponding  term  V*(H*u)  models  the  average 
conversion  of  viscous  work  of  the  fluid  into  two  quantities,  heat  and 
kinetic  energy. 

The  average  internal  energy  equation  is  derived  in  a  similar  fashion  as 
the  average  gas  continuity  equation  and  gas  momentum  equations.  We  multiply 
Eq.  (3.3)  by  3g,  integrate  over  the  averaging  volume  V(x)  and  use  formulas 
(2.9)  and  (2.21)  to  obtain 


at 


/VV  f  VVV  #  V  V  V  x  v 

$g  pe  dV  +  V*  j  $g  pue  dV  =  -  J  gn  *p(u~u  )e  dS 
V  V  S  sp  sp 
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-  /  $g  pV*u  dV 
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+  /  3g  *  ,  dV  -  V*  /  3g  Q  dV 
V  V 


(3.48) 


-  / 


g  Q*  n  dS 
®  sp 


We  define  the  average  specific  internal  energy  e  similar  to  the  average  gas 
velocity,  that  is,  as  the  quotient  of  the  average  internal  energy  density 
pe  and  the  average  mass  density  p  : 


e  =  J 


£e 


(t,x) 


p(t,x) 


(3.49) 


Using  Eqs.  (3.13)  and  (3.49),  Eq.  (3.48)  can  be  written,  after  some 
manipulation,  as 
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(ape)  +  V*  (apeu) 


a  pV  •  u  +  a$  ^  +  a4>  ^  -  V  < 


(aQ) 


*  SG 

+  P  VG6 


<d>  ~  ^  f  gQ*ncr«  dS 


sp 


+  V*  [apeu  -  aj  peu|  ]  +  dV  -  a4>  L  -  a$T]  (3.50) 


+  V*  [aQ  ~  I  BgQ  dV]  +  [apV«u  -  /  0gpV*u  dV] 


*  SG  *  1  r  ,pi 

-  p  —  e  <d>  -  —  /  ged  dS] 


VG 


where  and  Q  are  the  constitutive  models  for  the  average  dissipation 

function,  turbulent  dissipation  function,  and  the  average  heat  conduction, 
respectively.  A  The  average  energy  release  by  the  propellant  during  burning 
is  denoted  by  e(t,x).  The  term  otpeu-q  jveu  ,  which  is  -  (1/VG)  /ftgpe'u1  dV, 
is  analogous  to  that  in  Eq .  (3.26).  This  term  is  zero  if  either  e?  =  0,  or 

V  V  V 

u?  =  0,  i.e.  ,  if  e  or  u  is  a  function  of  time  only.  However,  in  turbulent 
flows,  this  term  can  be  significant.  A  model  of  tljig  term  as  the  gradient  of 
the  energy  variable  is  given  by  Chbeci  and  Smith  .  In  interior  ballistics 
the  term  is  probably  large,  because  for  moving  and  burning  grains  the 
extrema  of  ef  and  uf  are  likely  to  correlate.  We  denote  the  model  of  this 

terra  by  QT.  The  term  J  (gQ*  n  /SG)dS  represents  the  average  heat  flux  into 

p 

the  particle  from  the  gas  and  is  modeled  by  the  correlation  <e>.  The  models 

. 

for  4>  ^  and  4>  ,  Q  and  Q-j.,  and  e  and  <e>  are  discussed  in  Sections  4.7.3, 
4.7.4,  and  4.7.8,  respectively. 


We  now  can  rewrite  Eq.  (3.50)  as 


14 


T.  Cebeoi  and  A.  Smith ,  AndUjsis  of  Turbulend  Boundarii  Layers,  Academic 
Press,  New  York ,  1974. 
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—  (ape)  +  V*  (apeu)  =  Tip  V'u  +  at,  +  atm  ■  V*  (aQ)  -(a Q^,) 

0  t  L/  X 


* 

+  P 


SG  ~  /-K  _  SG  ,-s. 
VG  <d>  VG  <e> 


+  {^G  ^  ”  *L(t’x)  "  *T<t»x>]  dVl 

-  {V*i-  /  0g[Q(t,O  -  Q(t  ,x)  ]  dV}  -  { V*  [  (a  [peu]  -  apeu)  -  aQj} 

VLr  y 

-  /  g[Q(t,C)*nsp  -  <e>]  dS}  (3.51) 

s 

p 

+  {-^r  /  gg[p(t,x)V*u(t,x)  -  p(t,C)V*u(t,C)]  dV} 

VG  V 

+  {p*f§[e<d>-4  /  SSi  dS] }  , 

s 

p 

where  the  terms  enclosed  by  braces  are  error-type  terms. 

The  first  four  error  terms  depend  on  a  model  and  are  discussed  in  the 
appropriate  model  section  (see  Section  4).  The  remaining  two  terms  can  be 
written  by  following  similar  analyses  to  these  in  the  average  gas  momentum 
equation  derivation  as 

—  /  gg  [pV*u  -  pV* u]  dV  =  a (t ,x)  p(t,x)  [V*u(t,x)  -  V« u( t  (x) )  ] 

VG  V 

(3.52) 

+  a(t,x)  V*u(t,5(x))  [p(t,x)  -  p(t,£(x))] 


and 


SG  r  *  v 

VG  te  <d>  - 

1 

SG 

/  g  e  d  dS]  | 
s 

p 

„  1  *  SG  | 

<  p  — —  ■ 

A 

[  e  1 

<d>  -  d|  +  d  max 

1  ^  VG  1 

l  | 

i 

(3.53) 

e(t,x)  -  e(t,^(c1))  |} 
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where  5(x)  is  a  point  in  V  (V  is  assumed  connected)  and  CJ  is  a  point 

gas  gas  i  r 

on  the  surface  Sp^. 

The  error  represented  by  Eq.  (3.52)  consists  of  two  parts:  the  error 
made  by  using  the  divergence  of  the  average  velocity  for  the  divergence  of 

the  local  velocity,  and  the  error  made  by  using  the  average  pressure 

correlation  for  the  local  pressure.  If  both  p  and  V»u  were  functions  of 

time  only,  the  error  would  be  zero.  If  the  term  is  not  negligible,  then  a 

correlation  that  models  the  average  fluctuations  of  pV*u  from  pV.u  must  be 
included.  Most  often  the  term  is  neglected,  but  a  model  may  be  necessary  in 
some  turbulent  flows.  The  error  generated  by  replacing  the  surface  integral 

/■V 

of  ged/SG  with  the  product  of  correlations  e<d>  ,  Eq.  (3.53),  also  consists 
of  two  parts.  The  first  involves  the  approximation  of  d  by  <d>  which  is 
discussed  in  Section  3.2.1.  The  second  is  small  if  the  fluctuations  are 
small  of  the  local  internal  energy  from  the  specific  internal  energy  of  the 
gas  at  flame  temperature,  e.  In  practice,  both  errors  are  assumed  small. 

^  a  * 

If  not,  a  correlation  which  models  the  fluctuation  of  ed  from  e<d>  over  the 
surface  of  the  grains  must  be  included. 

3.2.4  Derivations  of  the  Surface  Average  Equations.  On  the  surface  of 
the  particles,  the  ^average  normal  regression  distance  d  and  the  average 
surface  temperature  T  can  be  defined  according  to  the  definition  (2.26), 

*  * 

where  d  and  T  denote  £he  local  values,  respectively.  For  a  spherical 
particle,  for  example,  d  is  the  local  difference  between  the  original  radius 
of  the  particle  ^and  its  current  radius.  According  to  Section  2.2.5,  the 
variables  d  and  T  satisfy  the  differential  Eq.  (2.41)  so  that  the  average 
regression  distance  equation  is 


*  * 


* 

ad 


u.Vd  =  <d>  +  /  g(|S.  -  <d>)  ds} 


+  {fe  /  (3  -  3)  (I  - » 


sp 


+  ib  I  (1  -  3)  .  H 


and  the  average  surface  temperature  equation  is 


(3.54) 
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* 

9_T 

3  t 


+ 


u*VT  =  <T>  +  {|g  /  g(§7"  <T>)  dS} 


+  /  (T  -  T)  {*  -  ucJ.VYg  dS} 


s  py 


+tio  j . 


(3.55) 


v  d  s  ^  # 

where  u  =  r— ,  <d>  is  the  correlation  for  the  regression  rate,  and  <T>  is 
Spot 

the  correlation  for  the  rate  of  change  of  grain  surface  tenperature. 


The  last  three  terras  in  each  of  the  Eqs.  (3.54)  and  (3.55)  are  error 
type  terms.  The  first  error  terms  in  Eqs.  (3.54)  and  (3.55)  are  the  surface 
averages  of  the  fluctuations  between  the  local  values  and  its  corresponding 
correlation  values  of  the  regression  rate  and  surface  tenperature, 
respectively.  The  regression  rate  term  is  discussed  in  Section  3.2.1  and 
similar  error  estimates  and  comments  can  be  made  concerning  the  surface 
tenperature  term.  The  remaining  error  terms  involve  fluctuations  from 

formally  defined  averages.  The  last  terms  in  Eqs.  (3.54)  and  (3.55)  involve 

v  V 

*  * 

fluctuations  of  d  and  T  from  their  average  values,  respectively.  Because 
the  integrands  of  these  surface  integrals  include  other  terms,  these 

integrals  are  not  surface  averages  of  fluctuations,  and,  thus,  are  not 

necessarily  zero.  The  other  set  of  error  terms  include  the  product  of  the 
fluctuations  of  the  local  interface  velocity  from  the  volume  average 

*  Ik  Ik 

particle  velocity  u  with  the  fluctuations  of  T  and  d  from  their  average 
values.  As  before,  the  integrals  involving  these  products  are  not  surface 
average  integrals.  If  the  fluctuations  are  small  over  the  surface  of  all 
the  particles,  then  the  terms  can  be  neglected.  Such  cases  occur  when  the 
regression  distance  and/or  the  surface  temperature  of  all  the  grains  are 

equal.  If  these  surface  integrals  represent  significant  contributions  to 
the  rate  of  change  of  the  variables,  correlations  for  them  must  be 

formulated  and  included  in  the  governing  Eqs.  (3.54)  and  (3.55). 


3.3  Summary  and  Discussion  of  the  Conservation  Equations  Without  Error 

Tterms 

In  this  section,  we  will  list  and  discuss  the  equations  derived  in 
Section  3.2  without  error  terms.  We  are  aware  that  some  of  the  neglected 
terms  may  be  significant  in  some  flows.  In  such  cases,  it  (they)  can  be 
appended  to  the  governing  equation(s)  and  modeled.  A  good  way  to  decide 
whether  a  term  should  be  neglected  or  included  in  a  set  of  equations  is  to 
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compare  the  accurate  solution  of  the  equations  with  data  from  well-defined, 
carefully  done  experiments.  Furthermore,  we  realize  that  some  of  the 
constitutive  laws  and  correlations  quite  possibly  can  be  coupled  to  each 
other  and  terms  in  the  governing  equations  could  be  grouped  differently. 
Thus,  the  formal  and  physical  meaning  of  some  of  the  constitutive  laws  and 
correlations  can  change.  Therefore,  the  form  of  the  equations,  correlates, 
and  constitutive  laws  for  interior  ballistic  applications  listed  in  this 
report  should  not  be  considered  as  final. 

The  porosity  Eq.  (3.20)  (the  average  solid  phase  continuity  equation) 
can  be  written  as 


(i-cO  +  v*  [(!-«)$]  =  -r  l 


(3.56) 


where  the  source  term  is  given  by 


r 


l 


SG 

VG 


<d> 


(3.57) 


The  average  solid  phase  momentum  Eq.  (3.44)  expresses  the  conservation  of 
the  solid  phase  momentum  density,  and  is 


9 _ 

9 1 


[ ( l-a)pu] 


+  V 


[  (  l-a)puuj 


-  (l-a)Vp  +  (l-a)p  A 


stress 


** 


(3.58) 


+  (l-a)p  Adrag  -  pur  j 


where 


( 1-a)P  Astress  =  +  (l-a)nT] 


(3.59) 


and 


(  l-«)p 


A, 

drag 


(3.60) 


The  average  gas  phase  continuity  Eq.  (3.21)  is 
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(3.61) 


|-£-  (ap  )  +  V*  (apu) 


The  average  gas  phase  momentum  Eq.  (3.30)  expresses  the  conservation  of  the 
momentum  density  and,  with  the  definition  of  drag  (3.40),  can  be  written  as 


|y  (apu)  +  V»  (apuu)  =  -  aVp  +  ap  Avigc  +  ap  \urh 


** 

+  PuT1 


(l-a)p 


A. 

drag 


9 


(3.62) 


where 


ap  A  .  =  V*  (all), 

vise 

and 


(3.63) 


ap  At  ,  =  V*  (an  )  .  (3.6  4) 

turb  T 

The  average  gas  phase  energy  Eq.  (3.51)  expresses  the  conservation  of  the 
gas  phase  internal  energy  density,  and  is 

(ape)  +  V*  (apeu)  *  -  apV*u  +  a$>  ^  +  a'F  ^  +  pel^  ,  (3.65) 

where 

4  i  =  4^  +  4,p  ,  (3.66) 

and 


1  =  -  V*  (aQ)  -  —  <*e>  -  V«  (aQT) 


(3.67) 


45 


Sect.  3.3 

The  term  $  ^  contains  all  the  models  for  the  heat  dissipation  functions  and 
the  term  V ^  contains  those  for  the  heat  conduction  within  the  gas  and  to  the 
particles,  and  the  turbulent  heat  flux.  The  average  specific  energy 

A 

released  by  the  burning  of  the  propellant  is  denoted  by  e(t,x). 

The  governing  equations  for  the  surface  average  regression  rate  (3.54) 
and  for  the  surface  average  surface  temperature  (3.55)  are 

* 

9  d  *  *  • 

—  +  u*V  d  =  <d>  ,  (3.68) 

and 

: k 

||+.u*VT«<T>  •  (3.69) 


Because  the  left-hand  sides  of  these  equations  represent  material 
derivatives,  one  can  interpret  the  equations  as  state  equations  for  the 
surface  material. 

SG  • 

The  source  term  is  modeled  by  which'  appears  in  every  volume 

averaged  equation.  Recalling  the  definition  of  the  source  term 


fd(t,x) 


V 

-4  /  8  1  dS' 


d  >  0 


(3.70) 


we  see  that  the  model  must  be  zero  when  no  particle  is  burning  within  the 

averaging  volume  at  point  (t,x)  (regression  rate  d  is  zero).  When  no 
particles  exist  within  the  averaging  volume  we  want  the  value  of  the  source 
term  to  be  zero.  This  is  reasonable  because  for  the  case  of  uniformly 
regressing  particles,  the  integral  in  Eq.  (3.70)  approaches  zero  as  the 
porosity  approaches  one.  Furthermore,  th£  value  of  the  model  must  be  always 
non-negative.  Comparing  Eq.  (3.61)  and  p  times  Eq.  (3.56),  we  see  that  the 

value  of  the  averge  mass  flux  per  volume  added  to  the  gas  phase  is  exactly 

that  being  taken  away  from  the  solid  phase  within  the  averaging  volume.  The 
average  balance  can  also  be  seen  in  the  momentum  equations  and  involves  the 
momentum  flux  model  pul1  The  drag  force  per  volume,  D/VG,  is  also  balanced 
on  the  average  in  the  momentum  Eqs.  (3.58)  and  (3.62).  We  note  that  the 
model  for  the  drag  force  D  should  be  zero  when  no  particles  exist  (a=l)  in 
the  flow  because  then  the  drag  force  Eq.  (3.40)  is  zero  (Sp  has  zero  surface 

area) .  Appropriate  types  of  average  stress  tensors  are  also  inglu^ed  in  the 

average  momentum  equations.  The  average  stress  tensors  II,  IIT,  II,  II T  in  Eqs. 
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(3.62)  and  (3.58)  are  defined  via  volume  averages  in  contrast  to  the  surface 
average  definitions  of  the  source  and  drag  terms.  These  average  stress 
tensors  are  weighted  by  the  appropriate^ phase  volume  fractions,  i.e.,  by 
a  for  H  and  II  ^  and  by  (1-ot)  for  II  and  H^.  Consequently,  the  contribution  of 
the  average  gas  stress  tensors  to  the  change  of  gas  momentum  consists  of  the 
two  terms  otV<(II  +11^)  and  Vot*(lI  +11^)  .  An  analogous  statement  holds  with 
respect  to  the  solid  phase.  The  internal  energy  of  the  gas  phase,  Eq. 
(3.65),  is  augmented  by  the  source  terra  pel^.  The  appropriately  weighted 
volume  averaged  heat  dissipation  functions  and  (the  contribution  from 
turbulence)  are  grouped  together.  The  average  work  done  by  the  gas  pressure 
is  denoted  by  -pV*u  and  is  weighted  by  the  porosity.  The  average  heat  flux 

SG  # 

between  the  gas  and  the  solid  is  represented  by  —  <e>.  The  correlation  <e> 

should  be  positive  when  the  temperature  of  the  gas  is  higher  than  that  of 
the  solid,  negative  in  the  opposite  case  and  zero  when  the  temperatures  are 
the  same  or  when  no  particles  exist  in  the  averaging  volume.  The  average 
heat  conduction  in  the  gas  is  modeled  by  V*(aQ).  The  turbulent  heat  flux 
vector  is  modeled  similarly  by  V*(aQ  )•  The  last  three  terms  are  grouped 
in  one  term  The  surface  averaged  equation  for  the  average  regression 

distance,  Eq.  (3.68),  has  a  non-negative  valued  *  right-hand  side  represented 

o 

by  the  correlation  <d>.  The  governing  equation  for  the  average  surface 
temperature,  Eq.  (3.69),  has  a  right-hand  side  that  usually  should  have  the 

a 

same  algebraic  sign  as  <e>. 

The  limiting  case  of  no  particles  within  a  region  is  of  particular 
interest  in  interior  ballistics  applications  because  such  regions  do  exist 
inside  a  gun  tube.  The  other  limiting  case  of  no  gas  does  not  exist  in  our 
applications  and,  thus,  is  of  no  practical  interest.  In  the  case  of  no 
particles  (a=l),  the  set  of  conservation  equations  greatly  simplify.  The 
source  terms  are  zero  and  the  drag  and  interface  heat  transfer  terms  are 
also  zero.  However,  It  is  important  to  notice  that,  first,  the  gas  phase 
equations  do  not  reduce  to  the  local  equations,  Eqs.  (3.1)  through  (3.3). 
The  simplified  set  (a=l)  differs  in  form  from  the  local  equations  because  it 
includes  the  turbulence  terras,  that  Is,  7«(II^),  and  V*(Q  ).  This  fact 

reminds  us  that  the  resulting  set  of  equations  is  still  a  set  of  average 
equations  for  a  finite  averaging  volume  V.  Secondly,  even  if  the  averages 
of  all  the  products  of  fluctuations  were  zero  (no  turbulence),  then  the  set 
of  equations  for  the  gas  flow  would  have  the  same  form  as  the  local 
equations,  but  the  solutions  would  not  be  the  same  in  general.  This  is  so 
because  the  quantities  p,  u,  and  e  are  averaged,  and  their  initial  and 
boundary  conditions  are  not  the  same  as  the  initial  and  boundary  conditions 
for  p,  u,  and  e  in  general.  Thirdly,  if  we  let  the  averaging  volume  go  to 
zero  in  the  simplified  set  (with  a  =  l),  the  turbulence  terras  would  be  zero 
because  the  fluctuations  are  averaged  over  the  averaging  volume  which  has 
zero  volume.  In  this  case  (a*l  and  V(x)*0)  the  averaged  equations  reduce  in 
form  to  the  local  equations  and  the  initial  and  boundary  conditions  should 
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reduce  to  the  local  conditions.  Thus,  the  solutuions  of  the  tvjo  setg  would 
be  identical.  Fourthly,  in  the  case  of  a-1,  the  equations  for  d  gnd  T.^Eqs. 
(3.68)  and  (3.69),  are  homogenous  (<d>*  =  <T>  =  0)  but  a  value  of  d  and  T  can 
be  computed  from  these  equations  if  u  is  defined.  Although  these  values 
would  be  physically  meaningless,  they  allow  the  solution  to  be  computed 
numerically  everywhere  without  tracing  the  internal  boundaries  of  gas  and 
mixture.  Because  these  internal  boundaries  cannot  be  predicted  ahead  of 
time  in  a  two-dimensional  flow  field,  this  provides  a  distinct  numerical 
advantage.  Fifthly,  the  average  solid  phase  momentum  equation  is 
identically  satisfied  when  a=l.  Thus,  the  components  of  the  vector  u  cannot 
be  determined  froig  Eq.  (3.58),  and  the  numerical  advantages  discussed  with 
respect  to  d  and  T  are  lost.  In  fact,  when  an  implicit  numerical  algorithm 
is  used  tg  solvg  Eqs.  £3.56)  through  (3.69)  directly  for  the  variables 
p,  a,  u,  u,  e,  d,  and  T,  it  can  be  shown  that  the  matrix  equation  which  must 
be  solved  for  a  new  time  level  of  values  is  singular  (the  rank  of  the  matrix 
is  deficient)  when  a=l.  To  avoid  this  situation,  we  can  algebraically 
manipulate  the  porosity  and  solid  phase  momentum  equations  into  a  non- 
conservative  form  when  3u/3t  has  a  coefficient  one.  Then  the  components  of 
u  can  be  defined  everywhere.  Another  advantage  of  solving  for  the  values  of 
u  ,  d,  and  T  directly  from  their  governing  partial  differential  equations 
when  a=l  is  that  their  values  should  be  continuous  at  a=l  if  the  equations 
approach  a  non-singular  form  at  a=l. 

In  Section  4.1,  4.2,  and  4.3,  we  discuss  better  forms  of  the  partial 
differential  equations  and  another  choice  of  dependent  variables  for 
numerical  treatment. 
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4.1  Basic  System  of  Governing  Equations 

A  system  of  conservation  equations  for  average  flow  properties  was 
derived  in  Section  3.  One  obtains  an  equivalent  system  of  differential 
equations  by  solving  the  conservation  Eqs.  (3.56)  through  (3.69)  for  the 
time  derivatives  of  the  dependent  variables.  Let  the  ensuing  system  be 
called  governing  equations  of  the  flow.  It  consists  of  the  following  set  of 
equations 


3p 

3t 


-  V*  (pu)  -£■[(!-  a  )V*  u  -  (u  -  *)•  V  ( 1  -  a)]  +  (P  q  P )  T2 


3e 

3t 


* 

=  -  u*Ve  -  £-V*u  +  (— - — )  —  +  —  ($>  +  )  , 

p  a  p  2  p  1  1 


3u  /  n\  1„  1/  P  r  l~a  A  ,  A  i  A 

—  =  -  (u.V)u  -  -Vp  "(u-u)-r2-T  Adrag  +  Avigc  +  Aturb 


* 
3  u 

3t 


*  .  *  1  i; 

-  (u*V)u  -  -rVp  +  V  A 


P 


P_ 

* 

P 


+  A 

drag  stress 


(4.1) 


3a 
3 1 


=  v*  ( ( i-o  )u)  +  r. 


* 

3d 

3t 


-  u* Vd  +  <d> 


* 

3T 

3 1 


-  *• VT  +  <T> 


The  system  is  closed  by  a  number  of  correlation  models  that  will  be 
discussed  in  detail  in  Section  4.7.  Presently,  we  merely  give  a  short 
exposition  of  the  corresponding  terms  in  Eq.  (4.1).  The  listed  arguments  of 
the  correlation  functions  are  only  representative,  indicating  the  most 
obvious  dependences.  The  actual  models  may  depend  on  fewer  or  on  more 
arguments.  Also,  all  models  depend  implicitly  or  explicitly  on  the 
averaging  volume  and  on  the  averaging  weight  function. 

The  equations  of  state  enter  the  system  in  form  of  a  relation  for  the 
pressure,  viz., 

p  =  p(p ,e)  ,  (Pa)  .  (4.2) 
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The  mass  source  due  to  the  phase  change  by  combustion  is  presented  by 

* 

s  (d) 

T  =  (1  -  a)  -p  <d>  ,  (1/s)  ,  (4.3) 

vp(d) 

SG  *  £  £  £ 

where  we  define  —  =  (l-a)Sp(d)/v  (d),  and  v  (d)  and  sp<d>  are  the  volume 

and  surface  correlations,  respectively ,  for  propellant  grains  with  the 

regression  distance  d.  The  quantity  <d>  represents  the  regression  rate 

correlation.  Generally  it  is  a  function  of  the  type 

<d>  =  <d>  (p ,  | u  -  u I  ,  9p/8t)  ,  (m/s)  .  (4.4) 


The  heat  dissipation  is  modeled  by  the  function 


$  1  =  $ j(u,T,a ,u,d) 


(W/m3) 


(4.5) 


where  T(p ,e)  is  provided  by  the  equation  of  state  correlation.  The  heat 
conduction  is  represented  by  the  function 


-  Ijd,  VT,  V.VT,  <T»  ,  (W/m3)  .  (4.6) 

The  last  argument  of  4f  ^  in  Eq.  (4,6)  is  the  rate  of  change  of  the  grain 
surface  temperature,  which  may  be  modeled,  e.g.,  by 


<T>  =  <T>  (T,T,|u-u|)  ,  (K/s)  .  (4.7) 

The  term  ^rag  represents  the  acceleration  due  to  the  drag  between  gas 
and  particles 


*  .  * 


^drag  Airag  U)>^,T)  ,  (m/s  ) 


(4.8) 


The  velocity  governing  equations  contain  three  more  acceleration 
terms.  They  are,  the  acceleration  by  the  laminar  viscosity 
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(4.9) 


the  acceleration  due  to  turbulence 


Aturb  =  ^tarb  (T,Vu,V*Vu, . . .)  ,  (ra/s2) 


(4.  10) 


and  the  acceleration  due  to  intergranular  stress  and  solid  phase  turbulence 


A  =  A 

stress  stress 


*  k 

(ot  ,d  ,7u , .  .  •) 


(m/s2) 


(4.11) 


The  system  of  governing  equations,  Eqs.  (4.1),  is  for  numerical 
solution  more  advantageous  than  the  system  of  conservation  Eqs.  (3.56) 
through  (3.69)  because  none  of  the  Eqs.  (4.1)  become  identically  satisfied 
as  a+  1.  This  permits  one  to  carry  out  the  calculations  throughout  the 
interior  of  the  gun  tube  without  tracking  the  boundaries  of  regions  with  a  = 
1. 


We  can  further  improve  the  equations  system  by  selecting  a  new  set  of 
dependent  variables.  The  choice  of  the  new  variables  and  the  corresponding 
new  system  of  governing  equations  are  described  in  Sections  4.2  and  4.3, 
respectively. 

4.2  Choice  of  Dependent  Variables 

4.2.1  Particle  Number  Function.  If  the  source  term  Y^  is  computed 
using  Eq.  (4.3),  then  one  can  expect  numerical  difficulties  as  Vp(d) 
approaches  zero.  Interpreting  the  equation  physically,  it  is  plausible  that 
Jl  -  a  ~  v  ,  so  that  Y  ^  vanishes  at  the  limit.  However,  because  a  and 
d  (and,  consequently,  v^(d))  are  separate  variables,  their  numerical  values 
wilL,  in  general,  approach  the  corresponding  limits  at  different  times  and 
locations.  In  a  computer  program,  the  situation  requires  special  safeguards 
to  prevent  overflow. 

The  special  programming  can  be  avoided  if  the  number  of  particles  is 
introduced  as  a  dependent  variable.  This  can  be  done  by  different 
approaches.  In  one  approach,  one  assumes  that  the  governing  equations,  Eqs. 
(4.1)  for  a  and  d,  and  thg  source  term  correlation  (4.3)  hold  exactly.  Then 
the  number  of  particles,  m(t,x),  can  be  introduced  by  a  formal  definition  in 
terras  of  already  defined  functions.  In  a  second  approach,  one  avoids  the 
use  of  the  correlation  (4.^3)  and  defines  m(t,x)  concurrently  with  the 
particle  volume  function  v  (d)  such  that  the  equation  for  a  in  the  equation 
system  (4.1)  is  satisfied  approximately.  Finally,  one  can  define  ra(x,t) 
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by  a  specific  Reasonable”  formula  and  then  seek  to  determine  a 
corresponding  function  vp(d)  such  that  the  equation  for  a  is  approximately 
satisfied.  Each  of  the  approaches  requires  some  approximations.  The  last 
approach  has  the  advantage  that  it  provides  guidelines  how  to  chose  the 
particle  volume  function  vp(d). 

*  * 
We  start  with  the  first  approach  and  define  m  in  terms  of  a  and  v  (d) 

as  in  Eq.  (2.45)  by 


m(t,x)  =  VG  (l-ct)/Vp(d) 


(4.12) 


The  two  governing  equations  for  a  and  d  in  Eqs.  (4.1)  are,  if  the  definition 
of  T 2  by  Eq.  (4.3)  is  used, 


9  (Ha) 

3 1 


V* ( (  l-a)u) 


and 


s  (d)  # 

(  1-a  )  •  P-  v  <d> 

vp(d) 


* 

3d  *  * 

£=■  =  -  u*  Vd  +  <d> 
3 1 


(4.13) 


Next,  we  express  a  in  terms  of  m  and  vp  using  Eq.  (4.12),  and  obtain 


a 


m 


wVd) 


(4.14) 


The  expression  (4.14)  is  substituted  into  the  first  Eq.  (4.13).  After 
simple  manipulations,  whereby  the  relation 

dv  ( d)  * 

- -  s  (d)  (4.  15) 

dd  p 


is  assumed,  one  obtains  from  the  system  (4.13)  the  new  system 


* 

3m 

3 1 


-  V' 


*  * 
(m  u) 


* 

3_d 

3  t 


-  u* Vd  +  <d>  . 


(4.16) 
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Thus,  one  can  replace  the  two  governing  Eqs.  (4.13)  by  the  two  Eqs. 
(4.16)  and  the  relation  (4.14).  If  m  is  used  instead  of  a  as  dependent 
variable,  then  the  source  term  T  ^  in  the  equation  system  (4.1)  is  calculated 
by 


2 


TO  k  » 

-  fc  sp(d)  <d> 


(4.17) 


instead  of  using  Eq .  (4.3).  The  expression  (4.17)  has  no  numerical 
singularities.  In  addition,  the  new  Eqs.  (4.16)  are  simpler  than  the 
previously  used  set  (4.13).  Physically  interpeted,  the  first  Eq.  (4.16) 
means  conservation  of  the  number  of  particles,  independently  of  their  size, 
whereas  the  second  equation  governs  the  average  size  of  the  particles, 
independently  of  their  number  in  the  averaging  volume. 

•k 

The  weak  point  of  the  ^described  formal  introduction  of  m(t^x)  (the 
first  approach)  is  that  m  and  the  governing  equation  for  m  contain 
inaccuracies  that  depend  on  the  quality  of  the  formula  (4.3)  for  the  source 
term  T  In  order  to  make  #  the  definition  of  m  ^dependent  ^  of  these 

inaccuracies,  one  can  define  m  concurrently  with  vp(d)  and  sp(d)  by  the 
relation  (4.12),  which  we  write  in  the  form 


m(t  ,x)  v  (d)  =  /  (H)  gdV 
P  V 


(4.18) 


the  Eq  •  (4.15),  and 

m(t,x)s  (d)  =  /  g  dS  *  SG  .  (4.19) 

S 

P 

The  Eqs.  (4.15),  (4.18),  and  (4.19)  are  consistent  in  the  sense  that  Eq . 

(4.19)  is  a  consequence  of  Eqs.  (4.15)  and  (4.18). 

The  exact  expression  for  the  source  term  is 


r2  =h  f  8  d  ds  =  ™  d 


VG 


Therefore,  if  Eq .  (4.19)  holds 


r2=WSp(d)  d 


(4.20) 


(4.21) 
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If  we  also  use  the  exact  average  value  d  instead  of  the  correlation  <d>  in 
the  governing  equation  for  d,  then  one  obtains  from  these  relations  and  from 
the  two  Eqs.  (4,13)  by  formal  manipulation  as  above 


* 
£  m 

9  t 


V-  (mu) 


* 

3_d 

3  t 


*  *  . 
u-Vd  -  d 


(4.22) 


Eqs.  (4.21)  and  (4.22)  are  derived  without  any  simplifying  approximations 
for  the  source  term.  When  the  equations  are  incorporated  into  the  equation 
system  (4.  1)  for  numerical  solution,  then  the  average  d  will,  of  course,  be 
replaced  by  the  corresponding  correlation  <d>. 

* 

The  weak  point  of  the  second  approach  is  that  the  two  functions  m  and 
Vp  with  the  desired  properties  do  not  exist  in  general,  and,  therefore,  one 
has  to  use  functions  that  satisfy  the  Eqs.  (4.15),  (4.18),  and  (4.19)  only 

approximately.  The  non-existence  can  be  seen,  e.g.,  by  considering  the 

ratio  Sp/Vp,  which  according  to  Eqs.  (4.18)  and  (4.19)  is  equal  to 

S  (d)/v  (d)  =  /  g  dS/  /  (1-6)  g  dV  .  (4.23) 

S  V 

P 

Jhe  right-hand  side  of  Eq.  (4.23)  obviously  depends  not  only  on  the  average 
d,  but  also  explicitly  on  t  and  x.  Even  in  the  special  case  where  all 

A  * 

particles  are  equal,  i.e.,  d  =  d  =  constant,  the  ratio  depends  on  the 
position  of  the  grains,  i.e.,  explicitly  on  t  and  x.  On  the  o^her  hand,  if 
g  is  a  constant,  then^Eq.  (4.23)  can  be,  indeed,  a  function  of  d  only,  and  a 
proper  function  Vp(d)  might  be  found.  (Actually,  g  can  be  only 

approximately  a  constant  in  order  to  insure  the  differentiability  of  the 

average  flow  variables,  see  Section  2.4). 

Because  the  Eqs.  (4.15),  (4.18),  and  (4.19)  cannot  l^e  satisfied 

exactly,  one  might  as  well  define,  as  a  third  approach,  m(x,t)  by  a 

reasonable  formula  and  then  seek  such  a  function  v  (d)  that  satisfies  the 

P 

abo^e  mentioned  equation^  approximately.  (The  other  possibility,  to  choose 
Vp(d)  and  then  define  m  by  Eq.  (4.18)  amounts  to  the  definition  by  Eq. 

(4.12).  The  corresponding  m  bas  undesirable  limit  properties  when  some 
grains  in  the  averaging  volume  are  reduced  by  combustion  to  zero.) 

* 

Either  of  the  following  two  formulas  define  functions  ra(t,x)  with 
reasonable  limit  properties  : 


54 


-Ilf  I  SCIS} 


Sect.  4.2.2 
(4.24) 


1=1  pi  S^V 

m  =  I  |f-  /  gdv}  .  (4.25) 

i=l  pi  V.S2V 

x 

In  these  equations,  m  is  the  number  of  grains  or  grain  parts  in  V,  spi  are 
the  surface  areas  of  the  grains,  are  their  surfaces,  vp^  are  the 

magnitudes  of  their  volumes,  and  V.^  are  their  volumes.  The  contribution  of 
a  grain  that  is  reduced  to  zero  volume  is  g(5^(t)  -  x)  ,  where  5i(t)  is  the 

location  of  the  grain.  When  all  grains  are  reduced  to  zero,  then  either  of 

the  formulas  produces 


m 


★  m 

m(t,x)  =  l  g{£  (t)  -  x) 
i=l 


(4.26) 


If  all  grains  have  the  same  finite  size,  then  the  formulas  reduce  to  Eqs. 
(4.18)  and  (4.1£),  respectively.  Finally,  if  g  is  constant  then  the 
contribution  to  m  of  each  grain  that  is  completely  inside  V  is  one,  and  the 
contribution  of  a  grain  partially  in  V  is  less  than  one,  in  accordance  with 
its  location.  Only  for  constant  g,  and  all  grains  located  inside  V,  the 
function  m  is  independent  of  d.  Therefore,  the  factorization  as  postulated 
by  Eqs.  (4.18)  and  (4.19)  can  be  best  approximated  if  the  weight  function  is 
constant  over  most  of  the  averaging  volume. 

If  m  is  defined  by  either  of*  the  Eqs.  (4.24)  or  (4.2£),  then  one  may 

select  the  volume  correlation  vp(d)  to  fit  the  choice  of  m  .  The  surface 

area  correlation  s  (d)  is  then  obtained  by  the  formula  (4.15).  The 

*  P 

selection  of  vp(d)  is  discussed  in  Section  4.7.9. 

4*2*2  Pressure  Logarithm  and  Entropy*  The  equation  system  (4.1) 
contains  two  thermodynamic  quantities  as  dependent  variables,  namely,  the 
density  p  and  the  specific  internal  energy  e.  One  can  replace  this  pair  of 
variables  by  a  different  pair  of  thermodynamic  quantities  and  replace  the 
first  two  equations  in  Eq.  (4.1)  by  corresponding  governing  equations  for 
the  new  pair.  The  variables  can  be  chosen  such  that  the  new  system  of 
equations  is  better  suited  for  numerical  treatment. 

First,  we  notice  that  up  to  six  equations  contain  the  gradient  of  the 
pressure.  The  handling  of  the  gradient  terms  can  be  simplified  considerably 
if  the  pressure  p  itself  is  chosen  as  a  dependent  variable  instead  of  p. 
The  replacement  reduces  the  total  number  of  terms  in  the  equation  system. 

Second,  one  may  replace  e  by  another  variably,  e.g. ,  by  the  specific 
entropy  s,  the  specific  enthalpy  h,  or  the  temperature  T.  These  choices  do 
not  simplify  the  equations.  The  number  of  terms  does  not  change  if  s  is 
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used  instead  of  e,  but  it  does  increase  if  h  is  used  instead  of  e.  Choosing 
T  as  a  dependent  variable,  one  obtains  the  most  complicated  equations. 

Based  on  these  considerations,  we  have  chosen  s  as  a  second 
thermodynamic  variable.  First,  it  does  not  complicate  the  equation 
system.  Second,  s  is  proportional  to  the  logarithm  of  the  temperature, 
whereas  e  is  proportional  to  the  temperature  itself.  Therefore,  if  the  flow 
contains  large  temperature  variations,  its  representation  in  terras  of  s  is 
much  smoother  and  more  amenable  to  numerical  differentiation.  (One  can 
expect  large  temperature  variations  in  certain  interior  ballistics 
problems . ) 

The  relation  between,  s,  p,  and  T  is  for  Noble-Abel  gases 


s  =  £n(T)  -  A^  Jtn(p) 


(4.27) 


with  constant  A^  and  A2 .  The  Eq.  (4.27)  suggests  that  q=J?,n(p)  would  be  an 
even  better  choice  than  p  as  the  other  thermodynamic  variable.  If  q  is  a 
function  of  p  only,  then  this  replacement  does  not  introduce  any  new 
complications  in  the  governing  equations.  Our  final  choice  of  thermodynamic 
variables  is,  therefore,  the  specific  entropy  s  [J/(kg*K)]  and  a  pressure 
logarithm  function  q,  which  we  define  as 

q  ( P )  =  q^JlnCp/pp  +  l]  ,  (Pa)  (4.28) 

with  constant  q-^  and  p-^. 

The  first  two  equations  in  the  system  of  governing  Eqs .  (4.1),  if 

expressed  in  terms  of  s  and  q,  are 

8s 

3 1 


9q 
3 1 

where 


=  -  u*7s  +  ^B  +  Hf  +  (4>  +  <J0 


(4.29) 


e  1  P 

=  -  u*'7q  -  — —  fv*u  +  ■=—  Bi  +  —  f e  -  e  -  e  Hi  F  -  —2-  ($  +  ill) 
pL  TJeL  sJ  p 
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B  =  —  f  (1-a)  V*  u  -  (u-u)*  V  (1  -a  )]  , 

a  L 


H  =  [  (e  +  p/p  )  -  (e  +  p/p  )] 


r  (d)  <d> 

a  p  2  a  p  VG  p 


(4.30) 


$  = 


—  $ 
Tt>  1 


f  = 


±-, 

Tt>  l 


and 


pq 

3p 

dq 

Ps 

3p (p , s) 

8  s 

 a  e(p ,s) 

dp 

eq 

3  P 

dq 

e 

3  e(p , s) 

s 

3  s 

the 

derivation 

of 

(4.31) 


A  9 


q 

In 


which  can  be  obtained  from  the  second  law  of  thermodynamics  (Rind). 

Eq.  (4.31),  dp/ dq  =  p/q^  by  Eq .  (4.28),  and  the  derivatives  of  the 

thermodynamic  functions  are  modeled  by  the  equation  of  state  correlation, 
described  in  Section  4.7.1. 


4,3  Final  System  of  Governing  Equations 

The  governing  Eqs.  (4.1)  can  be  expressed  as  follows  in  term  of  the  new 
set  of  variables  that  were  introduced  in  Section  4.2. 


l5F.  Hund,  Einfuhnmg_  in  die  Theovetische  PhysiJ^  Bd.  4 ,  "Theovie  dev  Wavme,  " 
p,  135  ff,  Bibliogvaphisch.es  Institute  Leipzig,  1950. 
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3q_ 

3t 


e  i  p 

u*  Vq  -  - —  (V*  u  +  B)  +  —  (e-e-e  H)r  -  —  ($+f) 
P„  T  e  s  p 

q  q  q 


a  u 

a  t 


* 

a  u 
aT 


k 

ajn 

aT 


P  n  *  1  — p, 

(U.V)  u  -  j-  Vq  -  (u-u)r  -  —  Adrag  +  Av±sc  +  Aturb> 


*  -  *  Pq  p 

Vq  +  ^  ^drag  +  ^stress 


=  ~  (u*V)  u 


* 

P 


*  £ 
V* (m  u) 


>(4.32) 


* 

3d 

a  t 


JL>  JL. 

u*  V  d  +  <d> 


* 
3  T 

3  t 


«JL»  JL. 

u*VT  +  <T> 


with 


*  jj-  [(1-a)  V*  u  -  (u-u)*V(  1-a)] 


a  =  1  -  v  (d)  m/VG 
P 


H  =  ^  [  (e  +  p/p)  -  (e  +  p/p)] 


r 


*  * 
J_  p_  m 

a  p  VG 


s  (d)  <d> 
P 


(4.33) 


The  partial  derivatives  pg,  pq,  eg ,  and  eq  are  defined  by  Eq.  (4.31).  The 
derivative  pq  =  dp/dq  is  equal  to  p/qj  if  q  is  the  pressure  logarithm 
defined  by  Eq.  (4.28). 

Models  of  the  various  correlation  terms  in  Eq.  (4.32)  are  discussed  in 
Section  4.7.  Their  physical  meaning  is  as  follows:  Y  represents  the  mass 
source  due  to  combustion,  $  represents  the  heat  dissipation,  V  contains  the 
heat  conduction  terms,  e(s,p),  T(s,p),  and  p(s,p)  are  thermodynamic  state 

functions,  Adrag  is  the  acceleration  due  to  drag,  Ay^sc  is  the  acceleration 
due  to  vicosity,  Atur^  is  the  acceleration  due  to  turbulence,  Agtress  is 
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the  acceleration  due  to  intergranular  stress  and  solid  phase  turbulence,  <d> 
is  the  regression  rate  correlation,  <T>  is  a  correlation  for  the  ^heat 
conduction  between  gas  and  particles,  e  is  e  at  flame  temperature,  sp(d)  is 
the  average  surface  area  of  a  single  grain,  and  vp(d)  is  the  average  volume 
of  a  single  grain.  The  variable  <T>  enters  also  the  first  two  Eqs.  (4.32) 
as  an  argument  of  the  term  ^ . 

The  correlations  are  defined  in  terms  of  volume  or  surface  averages. 
Therefore,  the  models  of  the  correlations  should  be  different  for  different 
averaging  volumes  and/or  different  weight  functions.  However,  because 
experimentally  determined  correlation  models  are  usually  reported  without 
reference  to  any  averaging,  their  relation  to  specific  averaging  procedures 
are  difficult  to  determine.  Therefore,  the  influence  of  their  relationship 
on  the  overall  accuracy  of  the  interior  ballistics  model  has  not  been 
established. 


4.4  Regions  of  Definition 

According  to  Section  2.3,  the  average  quantities  describing  gas 
properties  are  defined  at  all  interior  points  of  the  gun  tube,  except  for 
boundary  regions  the  shape  of  which  depends  on  the  averaging  volume.  The 
average  quantities  are  the  density  ap ,  the  energy  density  ape,  and  the 
momentum  density  vector  apu.  Consequently,  all  other  quantities  that  are 
defined  in  terms  of  these  quantities  are  defined  in  the  same  regions.  Such 
quantities  are,  e.g.,  e,  u,  s,  q,  T,  etc.  The  ^porosity  a  has  the  same 
region  of  definition.  The  grain  number  function  m  also  can  be  defined  in 
the  same  region  by  using  the  extension  m  =  0  if  the  averaging  volume 
contains  no  grains. 

Average  quantities  describing  grain  properties  are  defined  only  at 
reference  points  for  which  the  averaging  volume  conta^jis  grainy.  ^Therefgre, 
the  set  of  average  conservation  equations  for  (l-a)pu,  (l-a)p,  d,  and  T  is 
not  defined  in  regions  without  grains  (see  Section  3.3).  By  a  reformulation 
of  the  conservation  equations,  we  obtained  in  Section  4.3  an  equivalent  set 
of  governing  equations  (4.32).  This  set  has  no  singularities  at  a  =  1  and 
it  enables  one  to  calculate  nominal  grain  properties  at  all  interior  points 
where  the  gas  properties  are  defined.  Therefore,  one  can  extend  the 
definition  of  average  grain  properties  as  follows.  The  grain  properties  are 
defined  by  the  averaging  integrals  (see  Section  2.2),  if  the  averaging 
volume  contains  grains.  In  other  regions,  the  grain  properties  are  defined 
as  the  solution  of  Eqs.  (4.32).  In  ijtejior  ballistics  problems  this 
definition  amounts  to  an  interpolation  of  u,  d,  and  T  across  regions  without 
grains.  When  the  grains  have  been  reduced  to  zero  volume,  one  can  still 
calculate  their  motion,  which  now  corresponds  to  a  so-called  "dusty  gas" 
model.  In  such  a  gas,  the  dust  follows  the  gas  flow  according  to  a  drag 
law,  but  it  does  not  influence  the  gas  flow  itself.  Using  the  set  (4.32)  as 
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governing  equations  one  obtains  regions  of  "dusty  gas"  where  in  >  0  and  v  (d) 
■  0.  In  regions  ^with  m  -  0  and  v  (d)  >  0  the  equations  provide  an 
interpolation  of  u,  T  and  d  in  space  and*  time  between  regions  with  grains. 

In  the  boundary  regions  discussed  in  Section  2.3,  none  of  the  average 
quantities  are  defined  and,  consequently,  the  differential  Eqs.  (4.32)  have 
no  meaning  in  these  regions.  Strictly  speaking,  one  should  provide  boundary 
conditions  for  Eqs.  (4.32)  at  the  boundaries  Z/2  away  from  the  tube  walls 
and  1/2  or  £/3  away  from  the  breech  and  projectile,  if  the  average  volume  Is 
defined  as  a  sphere  (2.47)  or  cylinder  (2.49).  The  meaning  of  the  solution 
of  the  equations  in  the  boundary  regions  is  not  obvious  if  one  prescribes 
boundary  conditions  on  the  solid  boundaries  instead.  Section  4.6  contains  a 
discussion  of  the  boundary  condition  problem. 


4.5  Initial  Conditions 

Typical  local  initial  conditions  for  the  local  dependent  variables  in 
interior  ballistics  problems  are  constant  state  conditions  over  the  entire 
region.  Because  averaging  of  a  constant  produces  the  same  constant,  the 
intial  averages  in  most  cases  are  simply  equal  to  the  local  values. 

Deviation  from  a  constant  initial  state  typically  involves  either  a 
porosity that  is  not  uniform,  or  a  non-uniform  grain  size,  i.e.,  a  non- 
uniform  d.  In  these  cases,  one  cannot  use  the  local  values  of  m  and  5  as 
initial  values.  Instead,  the  initial  profiles  must  be  computed  by  averaging 
the  local  values,  whereby  the  same  averaging  volume  V  and  weight  function  g 
should  be  used  as  for  the  correlation  models  and  boundary  conditions. 

In  regions  where  intially  the  grain  number  m  is  zero  one  has  to 
extrapolate  or  interpolate  the  values  of  u,  d,  and  T  .  The  initial 

grain  velocity  is  normally  identically  zero  and  one  can  use  u  =  0  for  the 
extrapolation.  Likewise,  the  initial  grain  surface  temperature  is  usually 
constant,  and  the  same  constant  can  be  used  for  extrapolation.  The 
regression  distance  may  not  be  constant  if  different  sizes  of  grains  are 
located  in  different  regions.  In  such  cases,  one  has  to  use  a  common  sense 
extrapolation  that  produces  a  smooth  initial  surface  d  (0,x). 

In  the  boundary  regions,  "correct"  initial  values  cannot  be  specified 
for  reasons  explained  in  Section  4.4.  The  proper  choice  of  these  initial 
values  depends  on  the  method  of  treatment  of  the  boundary  regions.  However, 
one  can  assume  that  any  reasonable  treatment  will  produce  uniform  values,  if 
the  local  function  values  are  uniform.  Therefore,  one  may  specify  in  the 
boundary  regions  the  same  uniform  initial  values  as  in  the  interior 
region.  If  the  initial  conditions  are  not  uniform,  then  one  has  to  design 
such  an  extrapolation  of  the  averages  to  the  boundary  that  is  consistent 
with  the  treatment  of  boundary  conditions. 
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4.6  Boundary  Conditions 

A  theory  that  could  provide  guidelines  for  the  formulation  of  boundary 
conditions  for  averaged  equations  has  not  been  developed.  Therefore, 
interior  ballistics  calculations  usually  are  done  with  plausible  ad  hoc 
assumptions  about  boundary  values.  In  this  section,  we  shall  outline  the 
requirements  for  a  boundary  condition  theory  and  suggest  a  possible  approach 
to  the  formulation  of  such  a  theory.  Because  the  theory  has  not  been 
developed,  we  shall  also  discuss  ad  hoc  boundary  conditions. 

Discussing  boundary  conditions  for  averaged  differential  equations  in 
confined  volumes,  we  have  to  distinguish  between  two  boundaries.  For  the 
purpose  of  the  present  discussions,  we  call  them  the  outer  boundary  and  the 
inner  boundary,  respectively.  The  outer  boundary  consists  of  the  solid 
walls  of  the  volume.  In  interior  ballistics  the  solid  walls  are  the  tube 
walls,  the  breech,  and  the  base  of  the  projectile.  The  inner  bundary  is  the 
limit  of  validity  of  the  average  differential  equations.  As  discussed  in 
Sections  2.3  and  4.4,  the  inner  boundary  is  located  a  finite  distance  inward 
from  the  outer  boundary.  The  magnitude  of  the  distance  depends  on  the  size 
of  the  averaging  volume.  If  the  averaging  volume  is  a  sphere  with  the 
diameter  £ ,  then  the  inner  boundary  is  located  £  /2  away  from  the  tube 
walls,  breech,  and  projectile.  If  the  averaging  volume  is  the  cylinder 
described  in  Section  2.3,  then  the  inner  boundary  is  £  /2  away  from  the  tube 
walls  and  £  /3  away  from  the  breech  and  projectile  base.  Let  the  region 
between  the  outer  and  inner  boundaries  be  called  the  boundary  region,  and 
the  region  inside  the  inner  boundary  be  called  the  inner  region. 

Classical  theory  for  the  discussion  of  necessary  boundary  conditions, 
well-posedness ,  and  existence  can  be  only  applied  to  the  inner  boundary. 
Gough  (1974)  presents  some  of  the  discussion,  implicitly  assuming  that  the 
conditions  on  both  boundaries  are  identical.  The  assumption  is  permissible 
if  the  size  of  the  boundary  region  is  small  compared  to  the  size  of  salient 
structures  of  the  flow  field.  Because  the  size  of  the  boundary  region  must 
be  large  compared  to  the  size  of  propellant  grains  (see  Section  2.3),  it  is 
generally  not  small  compared  to,  e.g.,  the  gas  boundary  layer.  For  interior 
ballistics  flows,  therefore,  one  cannot  assume  that  boundary  conditions  on 
the  inner  and  outer  boundaries  are  identical. 

Ehysical  boundary  conditions,  such  as  u~uwa^]_>  are  only  given  for  the 
local  gas  phase  functions  on  the  outer  boundary.  The  only  physical  boundary 
condition  for  the  particles  is  that  no  single  particle  can  penetrate  the 
wall.  In  addition,  one  may  also  formulate  collision  conditions  for  single 
particles  impacting  on  the  wall,  i.e.,  on  the  outer  boundary. 

A  boundary  condition  theory  for  average  equations  has  to  bridge  the  gap 
between  the  outer  and  inner  boundaries.  It  should  provide  a  complete  set  of 
boundary  conditions  for  the  average  quantities  on  the  inner  boundary  in 
terms  of  the  local  physical  boundary  conditions  on  the  outer  boundary. 
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One  possible  approach  to  the  problem  is  by  construction  of  a 
continuation  of  the  solution  into  the  boundary  region.  If  such  a 
continuation  is  established,  then  one  has  reduced  the  problem  to  the 
formulation  of  boundary  conditions  on  the  outer  boundary  only.  The  simplest 
method  to  obtain  a  continuation  is  to  define  it  as  the  solution  of  the  same 
differential  equations  and  correlations  that  are  valid  in  the  inner 
region.  Then  one  needs  only  conditions  on  the  outer  boundary  and  disregards 
the  existence  of  the  inner  boundary.  This  is  the  usual  approach  in  two- 
phase  flow  calculations.  It  has  the  deficiency  the  one  has  no  guidelines 
how  to  formulate  the  boundary  conditions  for  the  continued  functions, 
because  they  are  neither  the  local  functions  nor  the  average  functions. 

A  more  promising  continuation  may  be  obtained  by  changing  the 
definition  of  the  averages  such  that  it  includes  the  boundary  region.  This 
requires  that  the  averaging  volume  V  has  a  shape  that  depends  on  the 
position  x  of  the  reference  point.  The  conservation  equations  of  Section  3 
are  derived  under  the  assumption  of  a  fixed  size  and  shape  of  V.  The 
averages  defined  for  a  variable  V  satisfy  a  different  set  of  differential 
equations.  The  continuation  into  the  boundary  region  could  be  computed  by 
solving  Eqs.  (4.32)  in  the  inner  region  and  the  new  set  in  the  boundary 
region,  and  by  matching  both  solutions  at  the  inner  boundary.  The  boundary 
conditions  on  the  outer  boundary  then  represent  conditions  for  averaged 
functions  and  can  be  modeled  accordingly. 

Because  a  theory  of  the  described  type  is  not  available,  we  now 

formulate  ad  hoc  conditions  that  may  be  used  for  the  differential  equation 
system  (4.32). 

The  local  boundary  conditions  for  the  gas  are:  U  =  uwall,  a  condition 
for  the  temperature  prescribing  either  T  =  Twg-Q  or  3T/3n  =  (3T/3n)  wa^^, 
where  n  is  the  normal  to  the  wall,  and  the  mass  conservation  equation.  In 
the  spirit  of  interpreting  the  solutions  of  the  differential  equations  as 
averages,  one  would  not  directly  use  these  conditions  as  boundary 
conditions.  Instead,  some  interpolation  is  needed  that  reflects  the 

averaging.  We  propose  the  following  approach. 

Let  1/ 2  be  the  distance  between  the  inner  and  outer  boundary  and  let 
e  be  an  estimate  of  the  thickness  of  the  gas  boundary  layer.  Let  <j>  be  a 
function  with  prescribed  local  boundary  value  and  n^  be  the  unit 

normal  to  the  inner  boundary,  pointing  outward  with  respect  to  the 
interior.  We  then  use  the  following  boundary  condition  on  the  outer 
boundary  for  gas  properties 

^outerb  “  f-Y  (^innerb  +  (^innerb*ni)  +  wall  ]/("2  e)  *  (4.34) 
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Because  l  is  larger  than  a  particle  diameter  (see  Section  2.3),  the  boundary 
value  on  the  outer  boundary,  when  computed  by  Eq.  (4.34),  will  approach  the 
local  boundary  value  only  if  the  particles  are  small  compared  to  the 
thickness  of  the  boundary  layer  (e»£/2).  This  may  be  the  case,  e.g.,  when 
the  flow  of  wear  reducing  additives  is  investigated.  If  the  particles  are 
large  compared  to  the  thickness  of  the  boundary  layer  (£/2>>e),  then  the 
outer  boundary  value  given  by  Eq.  (4.34)  approaches  an  extrapolated  value 
from  the  inner  boundary. 

Eq.  (4.34)  may  used  to  determine  the  boundary  values  of  u,  and  T  or 
9T/9n.  The  average  gas  continuity  equation  may  be  used  to  close  the  set  of 
boundary  conditions  for  gas  properties. 

The  formulation  of  a  boundary  condition  for  the  average  particle 
velocity  presents  a  dilemma.  On  one  hand,  the  condition  should  prevent  the 
particles  from  penetrating  the  wall.  On  the  other  hand,  the  average 
particle  velocity  at  the  outer  boundary  may  very  well  point  into  the  wall, 
merely  indicating  an  accumulation  of  particles  within  the  averaging 
volume.  As  an  ad  hoc  measure  we  disregard  the  second  possibility  and 
suggest  for  the  average  particle  velocity  at  the  outer  boundary  the 
following  formula.  Let  uDE  be  the  solution  obtained  from  the  differential 
equation  system  (4.32)  at  the  outer  boundary,  be  the  velocity  of  the 

wall,  and  the  nornia 1  to  the  wall  pointing  outward.  Then  the 

outer  boundary  value  of  u  is 


* 

u 


(4.35) 


satisfies  the  condition 


The  resulting  uouterb 


^uouterb  -  ^all^’Vall  **  ® 


(4.36) 


which  on  the  average  prevents  the  particles  from  flowing  through  the  wall 
and  permits  at  the  same  time  the  particles  to  leave  the  region  near  to  the 
wall  or  projectile. 

The  quantities  m,  d,  and  T  are  computed  by  solving  the  corresponding 
governing  equations  at  the  outer  boundary. 

4.7  Models  of  Correlation 

4.7.1  Equations  of  State.  For  the  derivation  of  the  average  equations 
in  Section  3,  we  used  the  averages  of  two  thermodynamic  quantities,  namely, 
the  density  p  and  the  specific  internal  energy  e.  The  conservation 
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equations  contain  two  other  thermo dynamic  quantities,  the  pressure  p  and  the 
temperature  T.  (The  latter  enters  the  heat  conduction  term  and  may  be  also 
used  in  other  correlations.)  They  were  assumed  to  be  related  to  e  and  p  by 
equations  of  state,  i.e.,  by 


p  =  p(p  ,e) 


and 


(4.37) 


T  =  T(p  ,  e) 


Generally,  one  uses  in  Eq .  (4.37)  the  same  functions  that  hold  locally. 

This  introduces  errors  in  several  terms  of  the  average  conservation 
equations. 

As  an  example,  let  us  consider  the  error  term  in  the  average  momentum 
equation.  Ihe  error  made  by  approximating  the  volume  average  of  the  local 
pressure  by  the  first  equation  in  (4.37)  is  from  Eq .  (3.31),  Section  3.2.2, 

%  -  ~  V[o  (t  ,x)  (p(t,£(x))  -  p(t,x))]  .  (4.38) 


As  discussed  in  Section  3.2.2,  to  minimize  the  error  by  a  proper  choice  of 
the  function  p,  we  need  to  minimize  the  errors  in  the  functional  values  as 
well  as  in  the  gradient  values.  However,  the  pressure  function  enters  the 
equation  system  in  various  places  and  in  different  combinations.  Therefore, 
the  use  of  the  local  equations  of  state  is  probably  as  good  an  approximation 
as  any.  Correspondingly ,  one  also  uses  the  local  equations  of  state  when 
the  entropy  s  is  introduced  as  a  dependent  variable. 

All  thermodynamic  variables  (temperature ,  pressure,  density,  energy 
entropy,  and  enthalpy)  are  completely  determined  in  terms  of  two  variables 
if  two  "equations  of  state”  are  provided  by  postulate  or  measurement.  Using 
the  two  given  equations,  all  other  relations  can  be  derived  from  the  laws  of 
thermodynamics,  which  provide  the  following  three  systems  of  differential 
equations  (Hund):  ^ 


3cv(p  ,  T) 

3p 


]_ 

2 

P 


a  p(p  ,t) 

2 

3T 


C  -c  -  T  f3p(p,T)12  f3p(p,T)1 

P  cv  2^3T  '  do 

P 


(4.39) 
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(cp  and  cv  are  the  specific  heats  (J/(kg»K)  for  constant  pressure  and  volume 
respectively)  , 


3e(p,T) 

3p 


3e(p  ,T) 
3  T 


=  _  J_  r  T  9  P-feolL  _p) 

9  1  9  x 


=  c 


(4.40) 


and 


3  s  (p  ,  T) 

3p 


1  3  p  (p  ,  T) 

2  3  T 

P 


3  s  (p  ,  T)  =  1 
3  T  T  v 


(4.41) 


An  equation  of  state  that  is  often  used  in  interior  ballistics  is  the 
Jfoble-Abel  equation 


p(p,T)  =  77  T  — 


M  1-f)  p 


(4.42) 


where  R  =  8.3143  J/(mol»K)  is  the  universal  gas  constant,  M  (kg/mol)  is  the 
molar  mass,  and  n  (m^/kg)  is  the  covolume.  From  Eqs.  (4.39)  and  (4.42)  one 
finds  that  for  a  Nsble-^bel  gas 


c„  = 


cv(T) 


and 


(4.43) 


=v<T>  +  5 


Therefore,  in  order  to  completely  specify  the  gas,  one  has  to  provide,  in 
addition  to  Eq .  (4.42),  a  temperature  function  cv(T).  Alternatively  one  can 
specify  instead  of  cy(T)  a  function  cp(T),  or  a  function  y(T)  that  gives  the 
ratio  Cp/ cy  =  y(T).  In  the  latter  case,  the  specific  heat  functions  are 
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cp(T) 


1  R 


Y (T) -1  M 


=  yCt)  R 

Y (T) -1  M 
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(4.44) 


We  assume  that  Y  (T)  is  constant,  and  obtain  with  this  assumption  the 
functions  e(p ,T)  and  s(p,T)  by  integration  of  Eqs.  (4.40)  and  (4.41).  After 
some  manipulations,  one  can  express  the  quantities  of  interest  in  terms  of  p 
and  s,  as  required  by  the  system  of  governing  Eqs.  (4.32).  The  results  are 
listed  below. 

integration  constant  for  the  entropy 


and  p^  are  reference  values  which  determine  the 


Y-i 


R 


T(p,s)  =  TR(Jy  Y  exp  s) 


K 


e  =  T 
Y~1  M 


J/kg 


(4. 45) 


f  R  T  'i  -1 

p  -  (m  p  +  "J 


kg/m3  , 


9T(p,s)  _  Y-l  T 
9p  Y  P 


9  T (p  , s )  _  Y-l  M  T 
9  s  Y  R 


(4. 46) 


9  e (p  ,  s )  _  Y~1  £ 
9  p  Y  P 


9  e(p  ,  s)  =  j_ 
9  s  y 


(4.47) 


and 


9p  (p  ,s) 

9  P 


(4.48) 
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The  square  of  the  sound  speed  is 


a2  =  y  £  _L_ 

T  p  l— np 


2/. 2 


mz/ s 


(4.49) 


The  specific  entropy,  expressed  in  terms  of  pressure  and  temperature,  is 

1-Y 


J/ (kg*  K) 


(4.50) 


4.7.2  Acceleration  by  Gaseous  Stresses.  The  governing  equation  for 
the  average  gas  velocity  in  the  equation  system  (4.32)  contains  the  terms 
^  and  Aj. ur ^ .  The  former  term  represents  the  acceleration  due  to  laminar 
viscosity.  The  latter  term  represents  the  acceleration  due  to  turbulence. 
A  simple  turbulence  model  is  a  Reynolds  stress  model  with  viscosity 
coefficients  depending,  e.g.,  on  temperature.  Then  the  forms  of  A^gj.  and 
Aturb  are  identical  (see  Eqs.  (B.9)  and  (B.34)).  We  restrict  our  discussion 
to  the  term  A^sc.  More  complicated  turbulence  models  are  possible  (see 
Gibeling  et  aljr  but  will  not  be  discussed  in  this  report. 

According  to  Section  3.3,  Eq.  (3.63),  the  viscous  acceleration  term  is 


A,  *  =  —  V  (an  ) 

'Vise  ap 


(4.51) 


where  II  models  the  gas  volume  average  of  the  local  viscous  stress  tensor  ^ 
ff.  The  loca^  tensor  is  given  in  terms  of  the  strain  rate  tensor  E  by 
(Tsien,  p.  13) 


V  V  ^  .V 

n  =  2  y  E  +  (X 


2 

3  y) 


trace  (E)  I 


(4.52) 


where  p  and  X  are  the  shear  viscosity  coefficient  and  the  bulk  viscosity 
coefficient,  respectively.  Both  are  assumed  to  be  functions  of 
temperature.  The  strain  rate  tensor  is  defined  by 

E  =i  (Vu  +  (Vu)T)  .  (4.53) 
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The  modeling  of  the  average  viscous  acceleration  term  involves  models 
of  the  average  viscosity  coefficients  and  a  model  of  the  average  strain  rate 
tensor  E. 

The  models  of  the  average  viscosity  coefficients  are  purely 
empirical.  A  convenient  set  of  formulas  is  the  following  generalization  of 
the  so-called  Sutherland  formula  : 


T1-5 

11  (W  ""o  +111^ 


and 


X1-5 

X(T)  =  Xo  +  x  i  xT+t 


Pa*  s 


Pa*  s 


(4.54) 


The  generalization  consists  of  the  addition  of  the  parameters  y  and  X  , 
thereby  including  in  the  model  the  constant  functions. 


The  average  strain  rate  tensor  E  is  usually  modeled  by  applying  the 
local  formula  (4.53)  to  the  average  velocities.  Then  II  is  obtained  by  using 
Eqs.  (4.52)  without  the  tildes  and  (4.54)  with  temperature  T(p,e)  calculated 
from  the  average  values  of  p  and  e.  The  approximation  error  is  Eq.  (3.32), 
Section  3.2.2,  divided  by  ap ,  i.e., 


St 


/  eg  [ff(u,p,e) 
V 


n (u,p ,e)] 


dV 


(4.55) 


V  N* 

The  error  part  that  comes  from  the  replacement  of  p  and  e  and  p  and  e  is 
probably  smaller  than  the  uncertainties  of  the  empirical  formulas  (4.54). 
However,  the  error  part  that  comes  from  the  use  of  the  average  velocity  in 
Eq.  (4.53),  can  be  large  because  the  formula  involves  derivatives  of  the 
velocity  and  in  a  viscous  two-phase  flow  the  local  derivatives  can  be  quite 
large.  The  integration  in  (4.55)  does  not  necessarily  cancel  out 
corresponding  large  local  undulations  of  the  integrand.  An  empirical 
correlation  based  on  careful  experiments  certainly  could  enhance  the 
usefulness  of  the  described  model  of  the  viscous  acceleration  terra. 

4,7,3  Heat  Dissipation.  All  the  heat  dissipation  terms  are  denoted  by 
and  they  enter  the  governing  Eqs.  (4.32)  for  the  specific  entropy  s  and 
for  the  pressure  logarithm  function  q.  According  to  Sections  3.2.3,  3.3, 
4.1,  4.2,  and  4.3  the  term  $  models 

ss  dV  '  <4-56) 
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where  the  local  heat  dissipation  function  $  is  given  by  (Tsien*®,  p.  15) 


4>  =  2  y  trace  (E^)  +  (X  -  y  )  (trace  E)^  ,  W/m'1 


(4.57) 


E  is  the  local  strain  rate  tensor  defined  by  Eq.  (4.53),  and  y  and  X  are  the 
local  shear  and  bulk  viscosity  coefficients,  respectively. 


Usually  $  is  defined  in  the  same  fashion  as  the  equations  of  state 
(Section  4.7.1),  i.e.,  by  calculating  a  $  with  the  same  formula  as  $,  but 
using  the  average  quantities  instead  of  the  local  quantities.  The  modeling 
of  the  viscosity  coefficients  is  discusse^in  Section  4.7.2.  In  Cartesian 
coordinates,  the  formula  is  (Tsien  ,  p.  15) 


_1_ 

PT 


*(E) 


1  1  9u-i2 

fc[i»  (sr  + 


8u,  2  0  8u.  2 

■)  +  (X  -  4w)  (air)  ] 


8  x 


(4.58) 


whereby  summation  over  i  and  j  is  assumed. 

Even  without  considering  turbulence,  Eq.  (4.58)  likely  underestimates 
the  value  of  the  expression  (4.56)  because  local  undulations  will  generally 
greatly  increase  the  value  of  the  integrand.  If  a  difference  exists  between 
the  average  velocities  of  the  phases,  then  local  velocity  gradients  entering 
Eq.  (4.56)  are  particularly  large,  but  are  neglected  in  Eq.  (4.58). 

In  order  to  estimate  the  effect  of  local  gradient  variations,  we 
compute  the  heat  dissipation  term  in  a  linear  flow  field  superimposed  by  an 
undulation.  Particularly,  we  assume  the  following  velocity  components  in 
Cartesian  coordinates  : 


»  A 

ul  =  U  +  ^  x  +  u(x,y,z)  , 


Un  =  u(x,y , z) 


(4.59) 


uo  =  u(x,y,z) 


where  the  undulating  part  is  given  by  the  formula 
u(x,y,z)  =  U  sin  x)  sin  y)  sin  z) 


(4.60) 
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The  local  heat  dissipation  for  this  flow  field  is  according  to  Eq.  (4.57) 


$ 


~  [-5-  y  (<f>  2  +  2  +  cj)  2  +  2<t>  2  +  2cJ)  2  +  2cJ)  2) 

p  T  L  2  T  xx  r  yy  Y  z  z  Y  xy  xz  yz 


XX 


+  (X  -|y)  {  (4>xx  +<j>yy  +<l>zz)2] 

where 

Kx  =  ^t~2  cos  (it  x)  sin  (r-  y)  sin  z)  +  2  fr 

♦yy  "  U  f1'  2  Sin  (l1  X)  COS  (f2-  y)  Sin  (f1-  Z) 

♦  zz  =  uf1  2  sin  (|2-  x)  sin  y)  cos  (^-  z) 

♦  xy  =  U  I1  sin  (x+y))  sin  z) 

♦xz  =  U  |2-  sin  (|2-  (x+z))  sin  y) 

♦yz  ■  u  I1  sin  (f1  x)  sin  (f5-  (y+z>) 


(4.61) 


(4.62) 


Next,  we  assume  that  the  averaging  volume  is  a  cube  with  side  lengths 
nL  and  that  the  weight  function  g  is  constant.  For  that  case,  the  integral 
(4.56)  yields 


$  = 


PT 


[(£)2cfu+x)  +(£)2 


TT2  (5|i+|x)] 


(4.63) 


The  first  term  in  the  brackets  in  Eq.  (4.63)  is  the  contribution  of  the 
linear  field  to  <& .  The  second  term  is  the  contribution  of  the  superposed 
undulations.  One  sees  that  for  Au/(nL)  ®  U/L  the  contribution  of  the 
undulations  is  about  40  times  larger  than  that  of  the  linear  flow  field. 
Interestingly,  the  contribution  of  the  undulations  does  not  depend  on  the 
number  of  periods  in  the  averaging  volume,  but  only  on  the  amplitude  U  and 
wave  length  L.  The  example  shows  that  the  usual  approximation  of  $  by  the 
formula  (4.58)  can  be  grossly  in  error. 
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A  model  of  the  contributions  of  undulations  in  two-phase  flow  due  to 
the  difference  between  u  and  u  can  be  derived  in  the  same  manner  as  Eq. 
(4.63).  To  simplify  the  formulas  let  us  chose  the  coordinate  system  such 
that  the  x-direction  coincides  with  the  direction  of  u-u.  Then  the  velocity 
undulations  may  be  approximated  by 

Uj  =  (u-u)  sin  (I1-  x)  sin  y)  sin  z) 


where  D  is  the  distance  between  the  centers  of  the  particles. 

let  ra  be  the  number  of  particles  in  the  averaging  volume.  We  associate 
each  maximum  of  the  function  u.  with  a  particle.  Then  there  are  four 
particles  in  the  elemental  volume  and  m  =  4V/D  .  Therefore, 


(4.6  4) 


D  =  (4V/m)1/3 


(4.65) 


The  contribution  of  the  undulations  (4.64)  to  the  dissipation  function  is 
one  third  of  the  contribution  of  the  undulations  (4.59)  in  all  velocity 
coordinates,  as  can  be  verified.  Therefore,  a  reasonable  model  for  the 
contribution  due  to  velocity  differences  is 

<$>  =  ^  (u-u)2(^)2/3  tt2(|p+|x)  ,  W*  (kg* K)  .  (4.66) 

In  a  computer  prograp  where  m  and  V  are  not  available,  one  can  use  in  Eq . 
(4.66)  the  quotient  m/(VG)  instead  of  m/V  without  changing  the  magnitude  of 
<$>.  The  correlation  (4.66)  probably  gives  only  the  order  of  magnitude  of 
the  contribution  due  to  velocity  differences  in  the  flow.  However,  it 
certainly  is  better  than  the  usual  assumption  <$>  =  0.  In  relation  to  the 
error  term  involving  the  dissipation  function  in  Eq.  (3.51),  Section  3.2.3, 
the  function  <$>  approximates  the  error  between  the  volume  average  of  the 
local  dissipation  function  and  the  average  dissipation  function  $(E). 

The  models  for  the  turbulent  dissipation  function  _vary  widely.  A 
simple  model  for  is  one  which  has  an  identical  form  £0  $  (Eq.  (4.58))  but 
with  different  viscosity  coefficients.  GLbeling  et  al.,  suggest  a 
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model  based  on  the  algebraic  relationship  among  a  turbulent  length  scale, 
turbulent  viscosity,  and  turbulent  kinetic  energy. 

The  complete  dissipation  term  that  enters  the  governing  equation  is  the 
sum  of  Eqs.  (4.58),  (4.66),  and  the  model  for 

=  Ip  *(E)  +  <$>  +  i=  *_  ,  W/(kg.K)  .  (4.67) 

pT  pl  i 


The  approximation  error  is  the  difference  between  the  expressions  (4.56)  and 
(4.67). 


4.7.4  Heat  Conduction.  The  heat  conduction  term  ¥  enters  the 
governing  equation,  Eq.  (4.32),  in  two  places.  The  term  itself  models  at 
least  two  phenomena:  the  heat  conduction  within  the  gas  defined  in  terms  of 
the  average  quantities,  and  the  heat  conduction  from  the  gas  to  the  solid. 

vw 

Depending  on  the  model  for  the  deviations  of  peu  from  peu,  we  also  can  have 
a  turbulent  heat  flux  vector  defined  in  a  similar  manner  as  the  average  heat 
conduction.  We  shall  discuss  each  of  these  models  in  turn. 

Locally,  the  heat  conduction  within  the  gas  is  assumed  to  be  governed 
by  Fourier's  law 


Q  =  -  k(T) VT  ,  W/m2 


(4.68) 


where  <(T)  is  the  thermal  conductivity  coefficient  which  depends  on  the 
local  temperature.  The  corresponding  average  heat  conduction  terra  in  Eq. 
(4.32)  is  a  model  of 


1 

apT 


/  6gQ  dv]  =  ^  V.[lg  J  eg<(T)  VT  dv] 


(4.69) 


The  volume  average  in  expression  (4.69)  is  usually  modeled  as  Eq.  (4.68) 

V 

without  the  tildes;  that  is,  the  local  temperature  T  is  replaced  by  the 
average  T  (obtained  from  the  average  values  of  s  and  q  by  the  equation  of 

state  correlation,  Section  4.7.1),  and  the  local  thermal  conductivity 

coefficient  k(T)  is  replaced  by  the  average  coefficient  k(T).  The  latter 

can  be  modeled  by  a  generalized  Sutherland-type  correlation, 


K(T) 


„1.5 


Ko  +  Kl  k„+T 


W/(m*  K) 


(4.70) 
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An  estimate  of  the  error  incurred  by  using  the  model  instead  of  expression 
(4.69)  can  be  obtained  as  follows  when  Vgag  is  connected: 


CQ  -  dv  -  aQ]  -  ^  /  Sg^T  dV  -  ck7T] 


I  A  A 

=  - —  V*  [a<V  T  -  aicV  Tl 

apT  L  J 


(4.71) 


A  A  A  A  A 

V  V 

where  T  =  T(s,q)  and  k  =  fc(T)  are  mean  value  points  of  the  integrand. 
Expanding  Eq .  (4.71)  further  one  obtains 


C 


Q 


1 

ap  T 


A  A 

V*  [a  (k-k)  VT  +  ukV(T-T)] 


(4.72) 


and 

IcJ  <  max  |— V*  [a  (k-k  )VT  +  a<V  (T-T)  I  .  (4.73) 

1  W'  'api  L  1 

v 

The  term  involving  the  difference  k-k  can  be  reduced  if  the  model  parameters 

k  ,  k  ,  ,  and  k  *  in  the  correlation  (4.70)  are  chosen  such  that 
o  1  l 

<(T)  ^  /  3g  k  dV  (4.7  4) 

The  term  involving  V (T-T)  reflects  the  modeling  error  due  to  local 
undulations  of  the  gas  temperature . 

The  heat  conduction  between  the  gas  and  the  particles  is  represented  in 
Eq .  (4.32)  by  a  model  of 


1  1_ 
apT  VG 


/  gQ*  n  dS 

S  sp 

P 


1  1 

apT  VG 


/  g  k  VT*  n  dS 

S  sp 

P 


(4.75) 


The  integrand  in  Eq .  (4.75)  is  the  heat  flux  into  the  particles.  We  define 
the  surface  averaged  heat  flux  by 


J-  / 

SG  1 


g  k  V  T*  n  dS 
sp 


W/  nr 


(4.76) 
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1  SG  • 
"  apT  VG  6 


(4.77) 


The  quantity  e  has  been  modeled  by 
correlations.  A  relatively  simple  formula 
* 

<e>  =  -2-  s  [h  (T-T)  +  h  (T-T)] 

SG  p  c  r 


various  different 
is  (Gibeling  et  al.) 


5 


W/m2 


experimental 


(4.78a) 


where  T  is  the  average  grain  surface  temperature.  The  coefficients  hc  and 
hr  in  Eq.  (4.78)  model  heat  transfer  by  conduction  and  radiation, 
respectively.  Gibeling  et  al.^  suggest  the  following  expression  for  the 
coefficients  in  case  of  spherical  particles  and  Nobel-Abel  gas  : 


h„  = 


* 

D  /2 
P 


+  0.2 


2  i  *  i  2 

r  y  R  (k2p)  |u-u|  >1 

T  y  D  /2 

P 


1/3 


* 

where  D  is  the  diameter  of  the  particles, 
coefficient  (Section  4.7.2),  and 


W/(m2K)  ,  (4.78b) 


and  p  is  the  shear  viscosity 


hr  =  *  crSB(TfT)  ( T2  +T2 )  ,  W/(m2K)  ,  (4.78c) 

where  e  is  the  particle  eraissivity  and  o  -  5.67032*  10  ®  is  the 

SB 

Stephan- Boltzmann  constant. 

.  . 

The  model  <e>  should  be  consistent  with  the  model  <T>  of  the  grain 

surface  temperature  rate  of  change.  The  relation  between  both  models  is 
discussed  in  Section  4.7.10. 


v  w 

The  model  of  the  significant  deviations  of  peu  from  peu  (denoted  by  Q^, 

see  Section  3.3)  can  have  different  forms.  One  model  of  the  turbulent  heat 

1 T  5 

flux  vector,  given  by  Ishii  and  Gibeling  et  al.,  is 

QT  =  -  kt[VT  -  |SL  (T±  -  T)]  ,  W/m2  ,  (4.79) 

where  T^  is  an  average  temperature  on  the  interface  (a  function  of  T  and  T) 
and  k^,  is  given  by  an  algebraic  formula  involving  an  effective  viscosity  and 
Prandtl  number. 


74 


Sect.  4.7.5 

The  heat  conduction  term  V  is  the  sum  of  the  three  described  models, 

i.e.  , 


^  ^  gas  +  ^particle  +  ^turb 

(4.80) 

■  Sr  7- <“VT)  -  f  <-e>  -  SSt  v- <aV  '  "/(k*'K)  • 


4.7.5  Acceleration  by  Drag.  The  acceleration  by  drag  between  gas  and 
particles  enters  the  governing  Eqs.  (4.32)  for  the  velocities  u  and  u.  The 
term  is  defined  by  (Section  3.3) 


Adrag  (l-a)p  VG 


(4.81) 


where  D  models 


W  I  g[nsp(P  -  P)  "  nSp*n]  dS  ,  (4.82) 

s 

p 

p  and  IT  are  the  local  pressure  and  viscous  stress  tensor,  and  p  is  the 

average  pressure.  In  interior  ballistics  applications,  the  term  is  modeled 
by  experimental  correlations  that  are  available # for  single  particles  (e.g., 
spheres)  and  for  packed  beds  of  particles.  For  situations  between  these 

extremes  one  has  to  interpolate. 

In  order  to  see  how  the  drag  coefficient  c^  for  a  single  sphere  relates 

to  A^rae;,  we  consider  a  situation  where  the  m  identical  particles  do  not 

interfere  with  each  other.  Then  the  absolute  value  of  the  drag  force  acting 
on  a  single  particle  is 


F 


i  1/  G-P) 

m  S 

P 


n  *11]  dS 
sp 


*  (l-a)  p  |A(jrag| 


(4.83) 


75 


In  terms  of  the  drag  coefficient  Cp, 


Sect.  4.7.5 

the  force  is  (Schlichting,  p.  15)^ 


F 


1 

2  CD 


*  1 2 

u-u  a  p 
1  P 


(4.84) 


where  ap  is  the  frontal  area  of  the  particle.  Eliminating  j F  J  between  Eqs. 
(4.83)  and  (4.84)  one  obtains 


I  *  I  1  |  *  1 2  m  1 

lAdragl  =  I  %  lU"Ul  apTCH  ’ 

or,  using  Eq.  (4.12),  Section  4.2.1, 
i  *  o  a„ ( ^ ) 

P 

The  drag  coefficient  for  a  single  sphere  can  be  approximated  by 


(4.85) 


(4.86) 


Cp  —  24/Re  +  0.4 


where 


Re  = 


*  I 
u-u 


D  (d)/p 
P 


(4.87) 


(4.88) 


* 

is  the  particle  Reynolds  number  and  Dp(d)  is  the 
diameter.  (About  the  approximation  (4.87),  see  Figure  1.5 

P.  16). 


average  particle 
in  Schlichting  , 


Substituting  the  expression  (4.87)  ijjito  (4.86)  and  observing  that  the 
acceleration  is  in  the  direction  of  u-u  one  obtains  for  non-interfering 
spheres  the  Reynolds  formula 


* 

(d) 


^Reynolds 


=  (u-u)  — 1 —  ( 0.2  | u-u |  + 


v  (d) 


12 


pDp(d)' 


(4.89) 


~^H~.  Rnhlichtina.  BnuodmiU  LOJd£l 1  Theory.  4th  Ed.,  McGraw-Hill  Book  Co., 
Hew  York,  1960. 
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e.g.,  the  Ergun  correlation  (Gibeling  et 


*Ergun 


*  ao(^ 

(u-u)  P - 


•k 

vp(d) 


4  \  (  1 . 7  5  I  u-u  I  +  150(l-a) 


3  2 
a 


P  Dp  (  2  ) 


(4.90) 


In  order  to  interpolate  between  both  formulas  one  may  assign  limits  for 
their  validity.  For  instance,  one  could  assume  that  the  dispersed  sphere 
formula  holds  for  a  >  0.9,  and  the  compacted  sphere  formula  holds  for  a  < 
0.65.  Then  the  acceleration  term  is 


A, 

drag 


4[(ot  -  0.65) 


^Reynolds 

^Reynolds 

^rgun 


for  a  >  0.9 

+  (0.9-a)A_,  ]  for  0.65<a  <0.9 

Dr  gun 

for  a  <  0.65 


(4.91) 


The  quoted  limits  are  arbitrary  and  may  be  changed,  if  experiments  are 
available.  Also,  other  than  the  Ergun  formula  may  be  used,  if  experimental 
data  indicate  a  better  approach. 

4.7.6  Acceleration  _Jz_  Granular  Stresses.  Acceleration  by  granular 
stresses  enters  the  governing  Eq.  (4.32)  for  the  particle  velocity  u.  The 
term  is  formally  defined  by  (see  Section  3.3) 

^stress - tT''[<1'«»‘l  ( - i-ipV.[U-a)i{T]  .  (4.92) 

(lnx)p  (l-a)p 


The  second  term  of  Eq.  (4.92)  represents  the  acceleration  of  the 
particulate  phase  by  solid  phase  turbulence  which  is  defined  by  Eq.  (3.46), 
£ection  3.2.2,  and  may  be  modeled  by  a  solid  phase  turbulent  stress  tensor 
II Because  the  density  of  the  solid  phase  is  much  larger  than  that  of  the 
gas  phase  and  the  sizes  of  the  propellant  grains  are  large,  the  turbulence 
of  the  solid  phase  is  assumed  negligible  and  IT  ^  is  set  equal  to  zero. 

* 

In  the  first  term  of  Eq.  (4.92)  the  variable  n  models 


1-ct  VG 


/  (1-0)  g  0T  +pl)  dV 
V 


(4.93) 
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(see  Eq .  (3,44)).  It  is  interpreted  physically  as  the  effect  of  grain 
interaction  with  grains  .  Without  such  an  interaction  the  stresses  II  inside 
the  grains  would  be  equal  to  the  negative  of  the  surrounding  gas  pressure 
or  nearly  so,  and  the  acceleration  term  Agtress  could  be  neglected,  except 
for  turbulence  considerations . 

Generally  in  interior  ballistics,  one  makes  two  assumptions  about  the 
model  II  of  the  average  integranular  stresses  .  First,  one  assumes  that  it  is 
a  function  of  a  only  and,  second,  one  assumes  that  it  is  a  diagonal  matrix 
i  .e  • , 


ff  -  I  R  (a) 
P 


(4  .94) 


The  second  assumption  means  that  the  stresses  have  the  effect  of  a  pressure 
that  acts  on  the  particles  in  addition  to  the  gas  pressure .  With  these 
assumptions,  one  obtains  from  Eq  .  (4.92)  for  the  acceleration 


^stress  i-a  da  [  *  Rp^a^  *  (4.95) 

P 

The  derivative  terra  in  Eq  •  (4.95)  is  interpreted  as  the  square  of  the  sound 
speed  a  in  the  dispersed  phase,  and  Agtresg  is  expressed  as 


^stress  a  2.  — oc  ^ 


(4.96) 


The  modeling  of  Astreae  is  reduced  by  these  assumptions  to  the  modeling  of  a 
sound  speed  function  a(a).  The  sound  speed  can  be  measured  in  packed  beds 
and  in  suspended  particle  flows,  so  that  the  model  can  be  tested  in  these 
special  cases . 

*  * 
The  function  a(a)  should  increase  with  higher  particle  density  (l-a)p, 

i  .e .,  with  decreasing  a.  Also,  as  a  approaches  one,  the  function  should 
approach  zero .  Let  agp  be  the  sound  speed  within  a  particle  and  ^ Let  us 
assume  th^t  for  a  =  a^  all  particles  touch  each  other,  so  that  a(a^)  = 
gSp  •  Let  a(a)  become  zero  at  a  =  1.  Then  a  reasonable  model  for 

a  (a)  is 


* 

a(a) 


a.-  a 

a  (— - -)  ( 

sp^a  -  a  J  K 
o 

0 


a0-a 

- — •) 

°2"°1 


for  a  <  a  <  a 
o 

for  <*2  <  a 


(4.97) 
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In  ^Eq •  (4.97),  the  value  a  -  a  corresponds  to  a  highest  density  ( 1- 

a  )p  that  can  be  achieved  by  compacting  the  particles.  If  a  *  0  then  one 
assumes  that  the  particles  can  be  crushed  and  compacted  to  a  solid  mass  with 
the  density  p.  The  last  factor  in  Eq.  (4.97)  merely  lets  a  approach  zero  as 
a  approaches  a^.  Thus,  one  assumes  that  for  a  >  particle  interaction  can 
be  neglected.  Gibeling  et  al.,“*  uses  a  similar  formula  in  which  «  0  and 
the  second  factor  is  set  equal  to  one.  Using  that  formula,  one  sets  a(a) 

=  0  for  a  >  o^.  It  seems  that  a  continuous  transition  to  zero,  as  provided 
by  our  formula  (4.97),  is  more  realistic. 


4.7.7  Burning  Rate.  The  burning  or  regression  rate  directly  enters 
the  governing  equation  for  the  regression  distance  d  in  Eqs.  (4.32).  The 
corresponding  term  is  defined  as  the  surface  average  of  the  local  regression 

*  * 

rate  3d/3t  =  (u  u)*  n  (see  Sections  3.1.2  and  (3.2.1)  and  is 

sp  sp 

approximated  by 

•k 

<d>  «  -|g  /  g  !y  ds  .  (4.98) 

s 

p 

The  linear  regression  rate  can  be  measured,  e.g. ,  in  closed  bomb  or 
strand  burner  experiments.  The  experiments  show  a  dependence  of  the  burning 
rate  on  the  gas  pressure,  on  gas  velocity  (erosive  burning)  and  on  the  time 
derivative  of  the  pressure  (dynamic  burning).  Best  established  is  the 
dependence  of  the  burning  rate  on  pressure,  which  is  modeled  by  the  equation 


B? 

d  =  B  +  B.p  (4.99) 

s  o  lr 

I  k  i 

with  constant  BQ ,  B^  ,  and  B2 .  The  dependence  on  the  relative  velocity  |  u-u  | 
and  on  the  pressure  change  3p/3t  can  be  incorporated  into  the  model  equation 

either  as  additive  terms  or  as  a  factor.  The  simplest  model  <d>  is  obtained 

•  • 

by  neglecting  these  dependence  and  setting  <d>  equal  to  d  ,  i.e., 

s 


for  <T>  <  T 

ignition 


<d> 


B„ 


J2  * 

B  +  B.p  for  <T>  >  T 
o  1  ignition 


(4.100) 


The  largest  uncertainty  of  this  model  comes  from  the  experimentally 
determined  model  parameters,  and  from  the  assumptions  that  erosive  and/or 
dynamic  burning  is,  or  is  not  important.  An  averaging  error  is  also 
introduced  by  the  use  of  the  equation  of  state  function  p(s,q)  in  Eq. 
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(4.100).  However,  the  error  is  likely  to  be  negligible  conpared  to  the 
general  inaccuracy  of  the  model  function.  These  errors  are  included  in  the 
error  estimate  (3.22),  Section  3.2.1. 

4.7.8  Source  Terms.  In  this  section,  we  discuss  terms  in  Eq .  (4.32) 
that  are  associated  with  the  burning  of  the  propellant.  They  are 
characterized  by  the  factor  <d>,  which  represents  the  regression  rate 
correlation  and  is  discussed  in  Section  4.7.7.  Because  of  this  factor,  the 
source  terms  are  equal  to  zero  if  no  burning  takes  place,  and  they  represent 
sources  of  mass,  energy,  and  momentum  if  the  grains  are  burning.  In  the 
governing  Eqs.  (4.32),  the  terms  have  the  common  factor  T  and  they  enter  the 
equations  for  s,  q  and  u.  The  factor  V  models  (Section  3.3) 


rSr^T:  J  ds 


a  p  VG 


sp 


sp 


and  is  defined  by 


r  «  ±  *-  •§£  <d> 

a  p  VG 


1/s 


(4.101) 


(4.102) 


In  Eq .  (4.102),  SG  can  be  eliminated  using  the  relation  (4.19),  Section 

4.2.1.  The  result  is 


r  * 


*  * 
i_  £L 

a  p 


vg 


<d> 


(4.103) 


as  stated  by  Eq.  (4.33). 

The  approximation  error  in  Eq .  (4.102)  is  that  of  the  correlation  <d> 
(see  Section  3.2.1).  In  the  expression  (4.103)  one  has,  in  addition,  errors 
associated  with  the  representation  of  the  weighted  surface  SG  the  product 
ms  .  Because  the  representation  is  part  of  the  definition  of  m  (see  Section 
4.2.1),  it  does  not  formally  introduce  new  errors. 

A  The  governing  Eq .  (4.32)  for  the  gas  velocity  contains  the  source  term 
(u-u)r .  The  term  models 

A  V  V 

^"^4?  /  (u_u)  g  [  (u  -u)*n  ]  dS  ,  rn/s2  .  (4.104) 

o  p  VG  g  sp  sp 

P 

The  error  in  the  governing  equation  caused  by  the  model  (4.103)  is 
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iP  1 
a  p  VG 


r  *  * 

J  (u-u)  g 

S 

P 


N/ 

* 

u)*  n 

sp 


].dS  • 


(4.  105) 


The  error  is  zero  if  all  grains  have  the  same  velocity  and  do  not  rotate. 

The  governing  Eq.  (4.32)  for  the  entropy  contains  the  source  term  HT . 
The  terra  is  derived  under  the  assumption  that  the  approximation 


.  V  V  tfc  A  -  V  tfc 

(  g  e  (u  -  u)*n  dS  ®  e  J  g  (u  -  u)»n  dS  (4.106) 

js  6  sp  sp  Js  6  sp  sp 

P  P 

holds.  Eq.  (4. 106)  is  indeed  an  identity  if  the  local  specific  energy  e  of 
£ he  gas  released  from  the  burning  propellant  surface  is  equal  to  a  constant 
e.  This  is  a  common  assumption  in  interior  ballistics.  The  constant  e  is 
the  specific  energy  of  the  gas  at  "flame  temperature",  i.e., 


y-l  M  flame  y-1  a  p 


g  I  ,  J/kg  , 


(4. 107) 


where  g  is  the  standard  acceleration  9.80665  m/s  and 
°a 


*p  Aflame  >  m 


(4.108) 


is  the  "force"  or  "impetus"  of  the  propellant.  (Sometimes  also  the  product 
gaIp(m^/s^)  is  called  the  "impetus"  of  the  propellant.) 

A 

In  some  cases,  a  modeling  of  e  may  be  better  than  the  assumption  of  a 
constant  e.  For  instance,  if  the  propellant  contains  a  retardant  then  one 
could  assume  that  the  flamg  temperature  is  a  function  of  the  regression 
distance,  and  consequently,  e  *  e(d).  Of  course,  the  modeling  then  involves 

A  * 

averaging  errors,  because  the  local  e(d)  would  be  replaced  by  a  function 
e(d)  of  the  average  d. 

The  factor  H  is  defined  by 


H  = -^  [  (e+p/p)  -  (e+p/p)]  ,  J/ (kg.  K) 


(4. 109) 


i.e.,  H  is  the  difference  between  the  enthalpy  of  the  gas  emerging  from  the 
flame  and  of  the  surrounding  gas,  divided  by  the  gas  temperature. 
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The  approximations  that  effect  this  term  are  those  of  the  equations  of  state 
(see  Section  4.7.1). 

The  source  term  in  the  governing  equation  (4.32)  Afor  the  pressure 
logarithm  function  q  has  as  a  factor  of  T  the  expression  (e-e-egH)/e  where 
eg(s,q)  and  eq(s,q)  are  the  partial  derivatives  of  the  specific  internal 
energy  with  respect  to  s  and  q,  respectively.  The  factor  is  derived  by 
formal  manipulation  and  the  approximations  involved  in  the  derivation  are 
the  same  as  discussed  above . 

4 .7 .9  Grain  Volume  and  Surface .  We  recall  the  discussions  in  Section 
4.2.1  about  the  definition  of  the  grain  number  function  m.  The  formal 
definition  of  the  average*  grain  volume  function  Vp(d)  and  of  the  average 
grain  surface  function  Sp(d)  should  be  consistent  with  the  definition  of 
m.  In  this  section,  we  shall  discuss  definitions  that  are  consistent  with 
Eqs .  (4.18)  and  (4.19),  respectively  . 

For  convenience,  we  repeat  the  pertinent  equations  and  definitions  in 
this  section.  Our  goal  is  to  find  such  functions  m,  v  ,  Sp,  that  satisfy 
the  following  set  of  relations 


(4  .110) 


(4  .111) 

m(t,x)  s  (d)  -  /  g  dS  (4  .112) 

1  S 

P 

m(t,x)  v  (d)  =/  (1-8)  g  dV  .  (4.113) 

p  V 

We  found  in  Section  4.2.1,  tha£  such  functions  in  general  do  not  exist  and, 
therefore,  suggested  to  define  m  by  either  of  the  following  two  equations: 


d(t ,x) 


dv  (d) 
P 

* 

dd 


SG 


g  d  dS 


* 

-  s  (d) 
P 


■  ■  l  /  g  dS} 

i=l  pi  s  nv 

p 

or 


* 

m 


m 


l  i 

i=l 


1 


Pi 


g  dv} 


(4  .114) 


(4  .115) 
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•k 

Once  ra  is  defined,  then  one  can  define  either  sp  or  Vp  by  Eqs.  (4.112)  or 
(4.113),  respectively,  and  find  the  other  function  from  Eq.  (4.111). 


The  approximations  involved  ^are,  first,  due  to  the  assumption  that  m, 
as  defined,  is  independent  of  d.  The  accuracy  of  the  approximation  is 
improved  if  the  weight  function  g  is  almost  constant  over  the  averaging 
volume.  A  second  approximation  is  due  to  the  assumption  that  Sp  or  Vp , 
defined  by  Eqs.  (4.112)  or  (4.113),  respectively  do  not  depend  explicitly  on 
t  and  x.  Again,  an  almost  constant  g  may  improve  the  accuracy  of  this 
approximation. 


The  modeling  of  the  functions  v^  and  sp  practically  is  done  at  a  limit, 
assuming  constant  g,  and  identical  particles.  In  this  case,  the  functions 
simply  represent  a  single  particle. 

If  there  is  a  variation  of  particle  sizes  within  the  averaging  volume, 
then  by  either  of  the  described  formalisms  one  obtains  average  n^  and  Sp 
that  are  slanted  towards  the  larger  particles.  Investigations  of  the 
significance  of  this  bias  have  not  been  done  for  interior  ballistics 
problems. 

4.7.10.  Grain  Surface  Heating  Rate.  The  grain  surface  heating  rate 
enters  the  governing  equation,  Eq.  (4.32),  for  the  grain  surface  temperature 

* 

3  T  *  *  . 

=  -  u* VT  +  <T>  .  (4.  116) 

o  t 


The  term  <T>  is  the  correlation  model  for 

f--§G/M*dS  >  (4-U7) 

P 

i.e.,  for  the  average  rate  of  change  of  the  surface  temperature.  The  change 
is  related  to  the  heat  flux  to  the  particles,  e,  discussed  in  Section 

4.7.4.  Therefore,  the  model  <T>  should  be  consistent  with  the  model  <e>. 

Like  the  grain  surface  and  grain  volume  functions,  the  surface 

temperature  model  function  is  usually  established  by  considering  the 
limiting  case  of  identical  grains,  i.e.,  by  treating  a  single  grain. 

Typically,  if  the  grain  has  a  simple  geometry,  one  calculates  the 
temperature  field  within  the  grain  corresponding  to  the  energy  transfer 
<e>.  This  involves  the  solution  of  a  differential  equation  that  is  valid 
within  the  grain  and  the  determination  of  the  corresponding  surface 
temperature.  However,  the  solution  of  Eq.  (4.116)  is  the  surface 

temperature  itself  and  the  temperature  field  within  the  grain  is  not  needed 
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(if  a  model  <T>  is  available.  Once  the  surface  temperature  is  known,  it  is 
used  to  determine  the  energy  transfer  at  the  next  time  step  or  the 
commencement  of  ignition.  This  type  of  calculation  is  used  if  one  is 
particularly  interested  in  the  ignition  process.  After  ignition,  all  heat 

transfer  is  assumed  to  be  zero,  because  then  the  energy  flow  phenomena  are 

dominated  by  the  combustion  and  the  associated  heat 
release.  The  continued  heating  of  the  grains  is  assumed  to  be  of  no 
consequence  for  the  combustion. 

In  order  to  illustrate  the  relation  between  <T>  and  the  heat  transfer 
from  the  gas  to  the  particle,  Section  4.7.4,  Eq.  (4.77),  we  consider  a  very 

simple  model  in  which  the  temperature  in  each  grain  is  assumed  to  be 

uniform.  (The  model  is  not  recommended  for  simulation  of  interior 
ballistics,  but  it  shows  the  salient  features  of  the  relation.)  Let  c  be 

p 

the  specific  £eat^  of  the  particle  material.  Then  the  heat  capacity  of  one 
particle  is  c  v  p,  (J/K),  and  the  relation  between  the  energy  transfer 
models  <e>  and  is 


m  0^* Vp<T>  -  <e>SG  .  (4.118) 


From  Eq.  (4.118)  and  expression  (4.77),  Section  4.7.4,  the  model  for  the 
heat  conduction  between  the  gas  and  particles  can  be  written  in  terms  of  <T> 
as 


* 

y  =  1  m 

particle  apx  VG 


[c  p*v  <T>] 
L  P  P  J 


(4.119) 


The  important  result  is  the  existence  of  a  relation  like  Eq.  (4.118)  between 
<T>  and  <e>.  It  would  be  replaced  by  a  different  relation  if  the  heat  flow 
within  the  particle  were  taken  into  account,  as  described  above.  In  that 
case,  the  expression  in  the  brackets  in  Eq .  (4.119)  would  be  changed 
correspondingly. 


5.  SUMMARY  AND  CONCLUSIONS 

Interior  ballistics  models  are  mostly  based  on  engineering 
approximations  and  insight,  like  Lagrange fs  model.  Alternatively,  one  can 
assume  that  the  gas  and  particles  locally  satisfy  all  conservation 
equations,  and  obtain  the  model  by  an  averaging  process.  In  this  report,  we 
follow  the  latter  approach  and  present  a  complete  mathematical  derivation  of 
weighted  volume  averaged  equations  including  all  error  terras,  sufficient 
conditions  for  the  necessary  differentiability  of  the  average  variables,  and 
regions  of  definition  of  the  average  variables.  Initial  and  boundary 
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conditions  that  are  consistent  with  the  volume  averaged  equations  are 
discussed.  Correlations  that  are  used  to  close  the  system  of  partial 
differential  equations  are  examined.  Some  of  these  correlations  are 
different  from  those  commonly  used  in  interior  ballistic  applications. 

The  average  governing  equations  that  are  derived  in  this  report  model 
the  transient  effects  of  viscosity,  heat  conduction,  and  turbulence  in  the 
compressible  gas  phase;  the  ignition,  intergranular  stress,  and  burning  in 
the  incompressible  solid  phase;  and  the  corresponding  interactions  between 
the  phases,  e.g.,  drag,  heat  transfer,  and  source  terms.  Turbulence  is 
defined  in  terms  of  volume  averages  with  only  elementary  models  presented 
for  completeness  of  the  report.  In  the  average  model,  quantities  appear 
that  are  defined  only  on  the  surfaces  of  the  grains.  We  show  that  these 
quantities  satisfy  a  general  partial  differential  equation.  The 

relationships  between  the  volume  average  equations  and  the  local  equations 

for  individual  phases  are  discussed  as  the  volume  approaches  zero.  Because 
these  equations  must  be  solved  via  a  computer,  an  appropriate  form  and 
choice  of  dependent  variables  for  numerical  solution  are  discussed.  Thus, 
this  report  presents  a  complete  and  consistent  mathematical  model  of 

interior  ballistics  for  non-reacting  burning  particle-gas  flows. 

The  exposition  of  the  theoretical  basis  of  averaged  equations  permits 

us  to  draw  the  following  conclusions: 

First,  the  proper  averaging  domain  is  a  finite  volume  that  is  larger 

than  the  propellant  grains.  Line  and  surface  averaging  cannot  be  used 

because  the  corresponding  averages  do  not  have  the  necessary 

differentiability  properties.  Infinite  volume  averaging  is  not  appropriate 

for  interior  ballistics  (or  other  confined  flows)  because  in  such  a  volume 
the  phases  do  not  occupy  complementary  spaces.  Time  averaging  is  not 

applicable  in  interior  ballistics  because  of  the  unsteady  and  rapidly 

changing  flow  conditions,  including  moving  boundaries. 

Second,  the  average  equations  are  valid  only  for  cases  where  the 

averaging  volume  consists  of  gas  and  particles  or  just  gas  and  where  the 

local  functions  have  no  discontinuities  within  their  respective  domains. 
Therefore,  average  governing  equations  are  not  suitable  for  describing  flows 
with  shocks,  contact  discontinuities,  etc.  On  the  other  hand,  by  a  proper 
formulation  of  the  governing  equations,  we  obtain  a  system  that  can  be 

solved  numerically  without  explicitly  following  the  boundaries  of  regions 

without  particles. 

Third,  the  average  equations  are  not  valid  in  the  boundary  region,  the 

thickness  of  which  is  equal  to  a  radius  of  the  averaging  volume. 
Consequently,  the  formulation  of  proper  boundary  conditions  is  problematic, 
and  has  not  been  solved  satisfactorily.  (We  do  not  consider  ad  hoc 

treatments  of  special  cases  an  adequate  solution  of  the  general  problem.) 
Also,  resolution  of  interior  ballistics  boundary  layers  based  on  volume 
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average  two-phase  equations  is  only  possible  in  exceptional  cases,  when 
particles  are  smaller  than  the  typical  boundary  layer.  The  study  of  two- 
phase  flow  fields  with  small  additive  particles  that  are  associated  with 
wear  and  erosion  of  gun  tubes  may  be  such  as  exceptional  case. 

s 

Fourth,  one-dimensional  interior  ballistics  models  based  on  volume 
averaging  are  less  problematic  than  two-dimensional  models,  because  the 
averaging  volume  occupies  a  finite  thickness  cross-section  of  the  tube  and 
is  large  compared  to  the  particles.  The  only  problems  with  such  models  are 
the  formulation  of  boundary  conditions  at  the  breech  and  projectile. 

Fifth,  a  mathematical  basis  for  two-dimensional  interior  ballistic 
models  could  possibly  be  obtained  by  an  extension  of  the  theory  of  average 
equations.  Such  an  extension  may  be  possible  by  generalizing  to  a  variable 
volume  averaging  or  by  using  statistical  averages  instead  of  volume 
averages.  The  first  approach  will  alleviate  some  problems,  but  it  cannot 
remove  the  basic  cause  of  problems  in  two-dimensional  modeling:  the 
particle  sizes  that  are  large  compared  to  the  gas  boundary  layer.  The 
second  approach  (statistical  averaging)  has  not  been  tried  successfully  for 
two-phase  flows.  There  the  encountered  problems  are  mathematical,  requiring 
a  major  investment  in  the  development  of  the  theory.  Also,  it  is  not 
certain  if  the  problems  associated  with  boundary  conditions  will  be 
alleviated  by  statistical  averaging. 
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APPENDIX  A 


App*  A 

GOVERNING  EQUATIONS  FOR  AXIALLY  SYMMETRIC  FLCWS  IN  CYLINDRICAL  COORDINATES 


This  appendix  contains  a  list  of  the  governing  equations  in  component 
form  in  cylindrical  coordinates  for  the  case  of  axially  symmetric  flow.  The 
subscripted  variables  denote  the  components  of  a  vector  and  not  the 
derivative  of  these  variables.  All  derivatives  are  written  in  a  non- 
abbreviated  form.  The  listed  equations  are  in  a  form  which  is  compatible 
with  Eq.  (4.32),  Sect.  4.3.  The  components  of  the  gas  average  velocity  and 
the  particle  average  velocity  are 


U  =  (ur,  Ug ,  Uz) 

,  (m/s) 

> 

(A.  1 ) 

^  ^  ^  ^ 
u  =  (ur,  V  uz  ) 

,  (m/s) 

(A.  2) 

where  the  subscripts  r,  9,  and  z  refer  to  the  radial,  angular,  and  axial 
coordinate  directions,  respectively.  The  components  of  the  gradient  of  a 
scalar  f  are 


7f  -  (ffU 


3fi 


The  divergence  of  a  vector  F  =  (Fr,  Ffl  ,  Fz)  is 


,  3(rF  )  3F 

u.F  =  I  _ _ £_  + _ £ 

r  3r  3  z 


(A. 3) 


(A. 4) 


The  independent  variables  are  time  t,  radial  position  r,  and  axial 
position  z.  The  dependent  average  variables  which  are  computed  from  the 
governing  partial  differential  equations  are:  the  specific  entropy  s,  the 

pressure  logarithm  function  q,  the  radial  gas  velocity  ur,  the 
circumferential  gas  velocity  ua ,  the  circumferential  particle  velocity  u. , 
the  axi^l  particle  velocity  u^,  the  number  of  particles  within  the  averaging 
volume  m,  ^the  regression  distance  d,  and  the  surface  temperature  of  the 
particles  T. 

The  entropy  equation  is 
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ii=-u  |^-u  p-  +  ^=B+Hr+$+f 
r  3r  z  3z  pT 


App.  A 

(A.5) 


where  p,  p,  T,  H,  T,  are  given  by  Eqs.  (B.6),  (B.  4),  (B.2),  (B.27),  and 

(B.26),  respectively.  The  expression  for  B  is 


&  A 

9  (ru  )  9u 


B 


If/,  \  ri  r  j_  z-\  ,  *  ^  9  (  1-a )  ,  *  9  ( 1-a  )i 

—  { ( 1-3 )  [T  —fr~  *  rrl  '  <ur  -  V  <v  V  S7p  • 

(A.6) 

and  the  porosity  a  is  given  by  Eq.  ( B. 1) •  The  dissipation  function  $  is 

$  =  $(E)  +  ^  $T  +  <*>  ,  (A. 7) 

where 


,  9u  2 

u  2  9u  2  u 

9  u  ; 

—  —  ii  r  r  + 

(  rl  +  f  z)  _  f  r 

r  + 

3  U  ^9r  J 

'  r  '  ^9z  '  *■  r 

9  r  I 

u„ 

_  9  u  9u  » 

9u_  0 

+  p  [ ( r  (— )) 

u-i)2 

l3r  3  r  1 

^3  z  ' 

3u  u 

3  U  0 

zl2 

9z  J  * 

(A.8) 


and  p,  X,  <$>,  and  are  given  by  Eqs.  (B.7),  (B.8),  (B.  13),  and  (B.35), 

respectively.  The  heat  conduction  term  ¥  is  given  by  Eq.  (B. 15)  as 


gas  particle  turb 


(A.  9) 


where 
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L 


i  r  i  3  (  3  t^  ,  3  (  aT>n 

gas  apT  Lr  9r  3  rJ  3  z  v 


3  zJ 


App.  A 

(A.  10) 


»turb  “  It  It  fi>  +lr(“,'ill 


■  1:  It  ( oct(ti_t)  It  -  h  (kt (TrT)  ffM 


(A.  11) 


3a  3 


3a  > 


Y  particle’  K  are  8iven  by  Ecls*  (B*17)>  (B*14),  respectively,  and  xr  Tt  are 
discussed  near  Eq.  (B.36). 

The  pressure  logarithm  function  equation  Is 


=  -  u  ia  -  u  ia_e_rl 


9(rur)  +  9uz  .  3e  1 


r  3  r  z  3  z  3p  A  r  3r  3z  3s  T 

3  q 


(A. 12) 


3  e 

3  q 


3s  3p_ 

3  q 


A 

where  p,  e,  T,  B,  e,  H,  T,  $  ,  and  'll  are  given  by  Eqs.  (B.  4),  (B.3),  (B.2), 
(A. 6) ,  (B.28),  ( B. 27) ,  (B.26),  (A.7),  and  (A.9),  respectively. 

The  radial  gas  velocity  equation  Is 

2 


3  u 


3  t 


=  -  u 


3  u  3  u  u.Q  j  i  ^  -  * 

— -  -  u  — -  +  JL  _  .Ip  i  -  (u-u  )  r 

r3r  z3z  r  dqp3r  rr 


(A.  13) 


(Adrag>r  +  <  AviJ  r  +  KurlPr 


where 


App .  A 


,  ~  0  3u  u  3u  ,  ~  3u 

(A  .  ]  =  —  {-r—  [ay  t  tt— - -  -  o— ■)  +  aX(—  (ru  )  +  3—^-)] 

v  visc; r  ap  l3r  L  3  v  3r  r  3z  '  vr  3r  r  3z  )j 


(A. 14) 


3  u_  3  u_  ^ 


3r^r 


and  p,  p,  T,  a,  y,  and  X  are  given  by  Eqs.  (B .  6 ) ,  (B.4),  (B.26),  (B.l), 
(B.7),  and  (B.8),  respectively.  The  radial  component  of  the  drag  )r 

is  given  by  the  radial  component  of  Eq.  (B.20). 
acceleration  due  to  turbulence  ^turb^r  cou^^ 


The  radial  component  of 
be  given  by  the  radial 


component  of  Eq.  (B.34)  which  is  Eq.  (A.  14)  with  y  and  X  replaced  by  y^  and 

V 


The  circumferential  gas  velocity  equation  is 


9ue  9ue  9ue  urue  ,  * 

Ft - ur  J7~  "  uz  3z“  ”  ~  ”  (u0"u0)r 


(A. 15) 


(1-a) 


(A,  )Q  +  (A  .  )Q  +  (At  ,  ). 


a  drag  0  vise  0  turb  0 


where 


^Avisc^0 


1_ 

ap 


{—  [ 

l3r  1 


ay  r 


3_ 
3  r 


(“)] 


+  h  fay  +  2lty  (■—)}  (A*16) 


3  u0 


3  z 


3_ 

3r 


and  a,  T ,  y,  X,  and  p  are  given  by  Eqs.  (B.l),  (B.26),  (B.7),  (B.8),  and 

(B.4),  respectively.  The  circumferential  component  of  the  drag  is 

given  by  the  circumferential  component  of  Eq.  (B.20).  The  circumferential 
component  of  the  acceleration  due  to  turbulence  (Atur^)Q  could  be  given  by 
the  circumferential  components  of  Eq.  (B.34)  which  is  Eq.  (A. 16)  with  y  and 
X  replaced  by  y^  and  XT. 

The  axial  gas  velocity  equation  is 


9uz  _  9uz  _  9uz  _  dp  1  3q  * 

3t  Ur  3r  Uz  3z  dq  p  3z  Uz  Uz 


(1-a) 


^Adrag^z  ^Avlsc^z  +  (A,-Iirt,) 


turb-'z 


(A. 17) 


where 
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.  .  3u  3u  3u  3u 

(Avisc)z  =  w  Iff  Mjr  +  rrh  +  ir  tar-  +  rH 


(A.  18) 


3  u 


3u 


8  u 


3  u 


♦<*(*r  +  jr£  +  rW 


and  p,  p,  T,  a,  y,  and  X  are  given  by  Eqs.  (B.6),  ( B.  4) ,  (B.26),  ( B.  1)  , 
(  B. 7)  ,  and  (B.8),  respectively.  The  axial  component  of  the  drag  (Adrag)z  is 
given  by  the  axial  component  of  Eq.  (B.20).  The  axial  component  of  the 
acceleration  due  to  turbulence  (^turb^z  could  be  Siven  by  the  axial 
component  of  Eq.  (B.34)  which  is  Eq.  (A. 18)  with  y  and  X  replaced  by  yT  and 

XT  # 

The  components  of  the  solid  phase  velocity  equation  are  the  radial 
solid  phase  velocity  equation 


* 

3  u 

i 

3  t 


3  u 


r  3r 


* 
3  u 


* 

ue2 


z  3  z 


-  dP  iLi  +  £_  (  a  )  +  (A  ) 

da  *  3  r  *  v  drag'r  ^stress'r  » 

4  P  P  (A.  19) 


the  circumferential  solid  phase  velocity  equation 


* 

9ue 

3  t 


* 

=  -  u 


3  un 


r  3  r 


* 

-  u 


* 
3  u, 


0 


z  3  z 


*  * 

UrU0  ,  p  (A  v 
r  +  *  ^Adrag^e 


(A. 20) 


and  the  axial  solid  phase  velocity  equation 


* 

3  u 
_ 2 

3  t 


* 
3  u 


* 
3  u 


=  -  u 


r  3r 


-  u 


z  3  z 


dP  _1  LSI  +  £.  (Aj  )  +  r a_  ) 

dq  *  3z  *  v  drag'z  ^^stress'z 


(A. 21) 


where  p  and  p  are  given^by  Eqs.  (B.6)  and  ( B.  4) ,  respectively.  The  density 
of  the  solid  phase  p  is  assumed  constant.  The  components  of  the 
accelerations  due  to  drag,  ,  and  intergranular  stress,  Agtress,  are 

given  by  the  components  of  Eqs.  (B.20)  and  (B. 23),  respectively. 


The  particle  number  equation  is 


* 
3  m 

Ft 


1  L_ 

r  3  r 


**  p) 

(rmu  )  -  — 
r  3  z 


** 

(mu  ) 
z 


(A. 22) 
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The  regression  distance  equation  is 

A 

*  *  * 

*  3d  *  8d  , 

-  m  -  XX  — r-  —  U  -  +  <d> 

3 1  r  3r  z  3z  * 

(A.23) 

where  the  burning  rate  correlation  <<!>  is  given  by  Eq.  (B. 25). 


The  surface  temperature  equation  is 

*  *  * 

3  T  *  3  T  *  3  T  •  v 

—  =  —  u  —  -  u  —  <T> 

3 1  r  3 r  z  3z 

(A. 24) 

where  the  correlation  <tT>  for  the  rate  of  change  of  grain 
temperature  is  discussed  in  Section  4.7.10. 

surface 

* 
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App.  B 


CORRELATION  MODEL  FORMULAS 

This  appendix  contains  a  list  of  correlation  model  formulas .  The 
formulas  are  discussed  in  detail  in  Section  4.7.  The  terms  listed  in  this 
appendix  are  in  a  form  compatible  with  Eq.  (4,32),  Section  4.3,  and  Appendix 

A. 

The  porosity  or  gas  volume  fraction  (Section  4.2.1)  is  given  by 
a  =  1  -  Vp( d)m/VG  .  (B . 1 ) 

The  equations  of  state  (Section  4.7.1)  are 


*(».«>-*£-)  <T-1)/Ye*p  (f^s) 

FR  1 

,  K 

(B.2) 

e  '  r-i  S  T 

.  J/kg 

(B.3) 

p  -  f + -O'1 

,  kg/m3  , 

(B.4) 

a2  =  yH-  * 

P  1— TIP 

,  m2/s2  , 

(B.5) 

where  R  =  8.3143  J/(mol»K)  Is  the  universal  gas  constant,  M  (kg/mol)  is 
molar  mass  and  n  (m3/kg)  is  the  covolume.  The  pressure  logarithm  functioi 
is  defined  by  (Section  4.2.2) 

1 

q  =  qjflnCp/pj)  +  l]  ,  Pa,  or 

p  =  Plexp  -  l) 

ql 

,  Pa  . 

(B.6) 

The  shear  viscosity  coefficient  p  and 
X  are  (Section  4.7.2) 

the  bulk  viscosity 

coefficient 

T1’5 

"  '%  +  uii,2  +  t  •  Pa’s  • 

(B.7) 

T1'5 

x  =  Xo  +xi  A0  +  T  *  Pa,s  • 

(B.8) 
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The  acceleration  by  viscosity  is  modeled  fyy  (Section  A. 7. 2) 

Ayigc  “  “  v  *  {a[2yE  +  (X  -  y  )  (trace  E)  i] }  ,  m/s2  ,  (B.9) 

where  E  is  the  strain  rate  tensor  computed  using  the  average  velocities, 
i.e. , 

E  -  0.5  (Vu  +  (Vu)T)  .  (B.  10) 


The  heat  dissipation  function  term  is  modeled  by  (Section  4.7.3) 


$  =  i_  $(E)  +  <$>  $T  ,  W/ (kg*  K) 


(B.  11) 


where 


$(E)  =  2y  trace(E^)  +  (X  -  -j  y  )  (trace  E)2  ,  W/m^ 

* 

v  1  |  *  1 2r  m  \2/3  2r  5  ,  1  .  \  T7/n  v\ 

<$>  =  _  |u-u|  ( 4  ;  vrJ  IT  (Ty  +TXJ  ,  w/(kg.  K) 


,  ( B. 12) 

,  (B.  13) 


and  is  given  by  Eq.  (B.35). 

The  thermal  conductivity  coefficient  k  is  modeled  by  (Section  4.7.4) 


t1.5 

<  -  K  +  K,  - ; — - 

o  1  ka  +  T 


,  W/(m-K) 


(B.  14) 


The  heat  conduction  term  in  the  governing  equations  is  modeled  by  (Section 
4.7.4) 


y  =  y  +  y  +  y 

gas  particle  turb 


where 


gas 


d? v-  (» * VI) 


W/ (kg*  K) 


(B.  15) 


( B. 16) 
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¥ 

particle 


! 


* 

_ 1_  m 

apT  VG 


s  [h 
pL  ' 


* 

(T-T) 


+  h  (T-T)] 


0 


before  ignition 
after  ignition 

(B.  17) 


with 


-T—+  0.2  (-Ijl  1/3  ,  W/(„2.  K) 


★  , 

D  /2 
P 


UDp/2 


and 


hr  =  *  aSB(»f)  (T2+T2) 


,  W/(n/*K) 


(B.  18) 


(  B. 19) 


^  __  o  _  o  _  A 

In  Eq.  (B.19),  e  is  the  particle  emissivity,  °gg  =  5.67032*  10  W*m  K  is 
the  Stephan- Bo ltzman  constant,  and  T  is  the  average  grain  surface 
temperature.  The  turbulent  heat  flux  within  the  gas,  is  given  by  Eqs. 

(  B.  36)  and  (B.37). 


The  acceleration  term  due  to  the  drag  between  gas  and  particles  is 
modeled  by  (Section  4.7.5). 


/  ^Ergun 

for 

a<0.65  , 

Adrag  =  |  ^a_0*65^AReynolds  +  (°,9-a  gun-1 

for 

0.65<a<0.9 

^Reynolds 

for 

0.9<a  , 

( B.20) 

where 

<1 

AErgun  =  («-«)  "T  t1*75  lU_*l  +  150  (1"a) 

P  « 

* 

,  m/s2  , 

(  B.  2  1) 
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^Reynolds 


-  (u-u)  [o.2  |u-u|  +  12  ,  m/s 

P 


* 

PD 


(B.22) 


The  acceleration  term  due  to  intergranular  stress  is  modeled  by  (Section 
4.7.6)  '  ~~ 


stress 


*2 

a 


1 

1-a 


V(l-a) 


m2/s2 


(  B.  23) 


where  a(a)  is  a  sound  speed  function  for  the  particulate  phase.  The 
function  is  modeled  by 


a  ,  -  a 


* 

a(a)  - 


(- — -)  ( 

^  n  —  of  J  ^ 


spva  -  a 


a0  -  a 

2 - ) 

a  2  al 


for  a  <  a  <  a 
o 

for  <  a 


(B .  24) 


The  burning  rate  is  modeled  by  (Section  4.7.7) 


<d>  -1 


0 


for  <T>  <  T 


I  B  +  B.pB2 
o  lr 


m/s ,  for  <T>  >  T 


ingition 

1 

ignition 


(  B.  25) 


The  source  term  T  is  (Section  4.7.8) 


r=HwsP<d>  *  • 


( B. 26) 


The  enthalpy  factor  H  of  the  source  term  (Section  4.7.8)  is  defined  by 


H  =  4  [  (e  +  p/p)  -  (e  +  p/p)] 


J/ (kg*  K) 


( B. 27) 


where  e  is 
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J/kg 


4 


1  R 
y-1  M 


Tf lame  y-l  g 


I 

a  p 
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( B. 28) 


with  g  =  9.80665  m/s2  being  the  standard  acceleration. 

* 

The  particle  geometry  enters  the  equations  as  the  four  functions  vp(d), 
s  (d),  D  (d)  and  a  (d).  We  provide  the  formulas  that  define  these  functions 
for  spherical,  cylindrical,  and  tubular  grains. 

* 

For  a  spherical  grain  with  initial  diameter  one  defines 


R  =  max  (0,  ((Dq  -  2d)/2) 


s  =  4  tt 
P 

ap  -  ,  fl2 


|  (B.29) 


A 

D  =  2R 
P 


* 

A  solid  cylindrical  grain  may  be  described  by  its  initial  diameter,  Dq, 
and  height ,  L *  Lat 


( B.30) 


L  =  L  —  2d 
o 


R  = 


*  * 

(DQ  -  2d)/2 


If  either  R  <  0  or  L  <  0,  then  the  grain  has  been  burnt.  If  both  quantities 
are  positive,  then  we  define 
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Vp  =  it  LR 


> 


sp  =  22ir  R(  R+L)  , 

ap  =  ( 2RL  +  it R2 ) / 2  ,  j  (  B.31) 

Dp  =  (2R+L)/2 

A  tubular  grain  may  be  defined  by  its  initial  height  and  the  initial 

outer  and  inner  diameters,  Dq  and  dQ,  respectively.  Let 

R  =  (Dq  -  2d)/2 

r  =  (dQ  +  2d)/2  ,  ^  B.32) 

*  * 

L  =L  -  2d 
o 

The  grain  is  completely  burnt  if  either  R-r<  0  or  L<  0.  If  both  of 

these  quantities  are  positive,  then  the  grain  geometry  functions  are 


VP  =  +  (R-r)L/2 


sp  =  *  (Dq  +  dQ)  (R-r+L) 
ap  =  (2RL  +  tt (R2  -  r2))/2 
Dp  =  ( 2  R+L) / 2  . 


( B.33) 


We  consider  a  detailed  study  of  turbulence  models  for  interior 
ballistics  flows  to  be  outside  the  scope  of  this  report.  Hence,  the 
correlation  models  are  quite  elementary  and  are  listed  in  this  report  only 
for  completeness.  The  acceleration  by  the  gas  phase  turbulent  stress  tensor 
^“turb  an<^  the  turbulent  heat  dissipation  function  could  have  the  same 

form  as  Ay^gc  (Eq.  (B.9))  and  $(E),  (Eq.  (B.12)),  respectively,  but 
different  viscosity  coefficients,  that  is, 


* 


* 
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,  (  B.  34) 


^urb  =  ^  V*  M  2  p  /  +  (X  T  "  t  p  J  ( trace  E)  XB  »  m/s2 
$T  =  2  yT  trace  (E2)  +  (XT  "  (trace  E) 2  ,  W/m3  ,  (B.35) 


where  \i  ^  and  XT  denote  the  viscosity  coefficients  for  turbulent  flows.  The 

manner  in  which  these  coefficients  are  determined  strongly  depends  on  the 

particular  turbulence  model  one  uses  and,  hence,  will  not  be  given.  As 

discussed  in  Section  4.7.6,  the  solid  phase  turbulent  stress  tensor  H_  is 

13  r 

set  to  zero.  The  turbulent  heat  flux  vector  QT  is  modeled  by  Ishii  and 
Gibeling  et  al.~*  as 

Qt  =  -  kt[VT  -  J2-  (T±  -  T)]  ,  W/m2  ,  (B.36) 

* 

where  T^  is  an  average  temperature  on  the  interface  (a  function  of  T  and  T) 
and  is  given  by  an  algebraic  formula  involving  an  effective  viscosity  and 
Prandtl  number.  The  corresponding  model  of  ¥turk  in  Eq.  (B. 15)  is 


^  turb 


1 

apT 


V*  (a  Qt) 


> 


W/  (kg*  K) 


( B.37) 


» 
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LIST  OF  SYMBOLS 


This  list  contains  symbols  that  are  frequently  used  in  the  report. 
Symbols  that  are  defined  and  used  only  locally  are  not  included  in  the  list. 

Function  symbols  in  general  indicate  average  quantities.  A  tilde  over 
a  function  symbol  is  used  to  indicate  the  local  value  of  a  function.  An 
asterisk  over  a  symbol  indicates  that  it  represents  a  property  of  the 


propellant  grains, 

► 

* 

a  - 

sound  speed  in  gas,  m/s 

asp 

* 

a  - 

sound  speed  of  particle  material,  m/s 

sound  speed  in  the  particulate  phase,  m/s 

aP 

2 

average  frontal  area  of  a  particle,  m 

adrag 

o 

acceleration  term  due  to  drag,  m/s^ 

^Er  gun 

2 

Ergun  correlation  for  A(jrag,*m/s 

^Reynolds 

2 

Reynolds  correlation  for  Acjrag,  m/s 

A  - 

stress 

2 

acceleration  term  due  to  intergranular  stress,  m/s 

cv 

specific  heat  capacity  at  constant  volume,  J/(kg*K) 

CP 

* 

d 

specific  heat  capacity  at  constant  pressure,  J/(kg*K) 

regression  distance,  m 

• 

d 

s 

stationary  burning  rate,  m/s 

<d> 

burning  rate  correlation  function,  m/s 

D 

drag  force  correlation,  N 

#  * 

D 

P 

average  particle  diameter,  m 

e  “ 

• 

specific  internal  energy,  J/kg 

A 

e  - 

e  at  flame  temperature,  J/kg 

es,eq 

partial  derivatives  of  e,  K  and  m'Vkg,  respectively 
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<e> 

correlation  for  surface  averaged  heat  flux  between  the 
particles  and  gas,  W/m^ 

E 

strain  rate  tensor,  1/s 

g 

averaging  weight  function 

H 

specific  enthalpy  difference  (4.109),  J/(kg*K) 

I 

» 

identity  tensor  of  second  order 

*P 

"force'*  or  impetus  of  the  propellant,  m 

l 

diameter  of  averaging  volume,  m 

m  - 

number  of  grains  in  averaging  volume 

* 

m 

weighted  number  of  grains  in  averaging  volume 

M 

molar  mass,  kg/mol 

n„  - 

sp 

unit  outward  normal  with  respect  to  the  gas  on  3p 

nsv 

unit  outward  normal  to  Sv 

p 

pressure,  Pa 

pq 

derivative  of  the  function  p(q) 

Q 

gas  phase  heat  conduction,  W/m^ 

Qt 

gas  phase  turbulent  heat  flux,  W/m^ 

q 

pressure  logarithm  function,  Eq.  (4.28),  Sect.  4.2.2,  Pa 

r  - 

radial  coordinate,  m 

R 

universal  gas  constant,  8.3143  J/(raol*K) 

% 

s 

specific  entropy,  J/(kg*K) 

sp 

average  surface  area  of  a  single  grain,  m^  * 

sp 

union  of  all  grain  surfaces  in  V 

SG 

weighted  area  of  Sp, 

Sv 

surface  of  averaging  volume  V 
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t 

time,  s 

T 

gas  temperature,  K 

T  — 

f  lame 

flame  temperature,  K 

* 

T 

grain  surface  temperature,  K 

<T> 

ft 

correlation  for  rate  of  change  of  grain  surface 
temperature,  K/s 

u  - 

gas  velocity,  m/s 

ur>u0>uz 

the  radial,  circumferential,  and  axial  components  of 
u,  m/s 

* 

u  - 

particle  velocity,  m/s 

•k  k  * 

VVuz 

£he  radial,  circumferential,  and  axial  components  of, 
u ,  m/  s 

usp 

velocity  of  a  point  of  S^,  m/s 

VP 

3 

average  value  of  the  volume  of  a  single  particle,  m 

V 

averaging  volume 

VG 

3 

weighted  value  of  V,  m 

x  - 

spacial  coordinate  vector,  m 

Z  -  ,] 

axial  coordinate,  m 

2 

surface  element  metric 

a  - 

gas  volume  fraction  (porosity) 

0 

phasic  function  (Section  2) 

* 

Y 

ratio  of  specific  heats 

.  r 

source  term,  (4.33),  Section  4.3,  1/s 

ri 

SG<d>/VG,  1/s 

r  2 

i/s 

P 

C 

surface  coordinate  vector 
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a* 


n 

O 

-  covolume  in  equation  of  state,  m  /kg 

< 

-  thermal  conductivity  coefficient,  W/ (m*  K) 

X 

bulk  viscosity  coefficient,  Pa*s 

y 

shear  viscosity  coefficient,  Pa* s 

5 

-  spacial  coordinate  vector,  m 

n 

-  viscous  stress  tensor,  Pa 

11  T 

* 

n 

* 

-  gas  phase  turbulent  stress  tensor,  Pa 

-  intergranular  stress  tensor.  Pa 

* 

11  T 

solid  phase  turbulent  stress  tensor,  Pa 

P 

o 

gas  density,  kg/m 

* 

P 

o 

particle  density,  kg/nr 

P  9  P 
s’  Kq 

-  partial  derivatives  of  p(s,q),  (kg/mJ) • (kg*  K/J) ,  and 

o  o 

s  /m  ,  respectively 

* 

-  function  describing  a  gas  property,  Section  2 

* 

<J> 

-  function  describing  a  particle  property,  Section  2 

$ 

dissipation  terra,  W/(kg*K) 

$  =  $ 

Li 

dissipation  function,  W/m^ 

$ 

T 

3 

-  gas  phase  turbulent  dissipation  function,  W/m 

<$> 

-  dissipation  correlation  term,  W/ (kg*  K) 

3 

-  p  T  $ ,  W/m 

\Kt,x,5) 

-  general  function.  Section  2 

heat  conduction  term,  W/ (kg*  K)  • 

gas 

-  heat  conduction  due  to  gas  conductivity,  W/ (kg*  K) 

ur 

particle 

-  heat  conduction  due  to  heat  loss  to  particles,  W/(kg*K) 

*turb 

turbulent  heat  conduction,  w/(kg*K) 
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data,  procedure,  source  of  ideas,  etc.) _ 


5.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far 


as  man-hours 
etc?  If  so. 

or  dollars  saved,  operating  costs  avoided  or  efficiencies  achieved 
please  elaborate. 

6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future 
reports?  (Indicate  changes  to  organization,  technical  content,  format,  etc.) 

Name 


CURRENT 

ADDRESS 


Organization 


Address 


City,  State,  Zip 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the 
New  or  Correct  Address  in  Block  6  above  and  the  Old  or  Incorrect  address  below. 


Name 

Organization 

ADDRESS 

Address 


City,  State,  Zip 


(Remove  this  sheet  along  the  perforation,  fold  as  indicated,  staple  or  tape 
closed,  and  mail.) 


